Activated and renewed HBC clustering and comparison with scATAC-seq

here::set_here()

## File .here already exists in /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/scripts

suppressPackageStartupMessages({
  library(slingshot)
  library(SingleCellExperiment)
  library(cowplot)
  library(rgl)
  library(clusterExperiment)
  library(RColorBrewer)
  library(ggplot2)
  library(pheatmap)
  library(wesanderson)
  library(gridExtra)
  library(irlba)
  library(tidyverse)
})

Import data

sds <- readRDS("../data/finalTrajectory/sling.rds")
counts <- readRDS("../data/finalTrajectory/counts_noResp_noMV.rds")
# counts <- round(counts)
sce <- readRDS("../data/finalTrajectory/sce_tradeSeq20200904.rds")
load("../data/ALL_TF.Rda")
clDatta <- readRDS("../data/finalTrajectory/dattaCl_noResp_noMV.rds")
timePoint <- readRDS("../data/finalTrajectory/timePoint_noResp_noMV.rds")

Select activated and renewed HBC clusters

I will select

• activated HBCs as cells sampled in the 24h timepoint that are in the starting cluster
• renewed (‘resting’) HBCs as cells sampled in the 14D timepoint that are in the HBC cluster

cols <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33", "#A65628", "#F781BF", "#999999")
view3d( theta = 10, phi = 10)
# Datta labels
plot3d(reducedDim(sds), aspect = 'iso', col=cols[clDatta], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")

# timepoints
plot3d(reducedDim(sds), aspect = 'iso', col=cols[timePoint], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")

# cluster labels used for TI
clSds <- apply(slingClusterLabels(sds), 1, which.max)
plot3d(reducedDim(sds), aspect = 'iso', col=cols[clSds], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")


actHBCId <- which(clSds == 8 & timePoint == "24HPI")
restHBCId <- which(clDatta == "HBC" & timePoint == "14DPI")
grpHBC <- rep(NA, length(c(actHBCId, restHBCId)))
grpHBC[1:length(actHBCId)] <- "activated"
grpHBC[is.na(grpHBC)] <- "resting"
grpHBC <- factor(grpHBC)
countHBC <- counts[, c(actHBCId, restHBCId)]

Dimensionality reduction

library(scran)
library(uwot)
library(scater)

## 
## Attaching package: 'scater'

## The following object is masked from 'package:clusterExperiment':
## 
##     plotHeatmap

logFracs <- log(sweep(countHBC+1e-5, 2, colSums(countHBC), "/"))
gv <- modelGeneVar(countHBC)
hvg <- getTopHVGs(gv, n=1000)
pc10 <- irlba::irlba(logFracs[hvg,], nv=10)
um10 <- uwot::umap(pc10$v %*% diag(pc10$d))

plot(um10, col=alpha(c("orange", "darkseagreen3")[as.numeric(grpHBC)], .75), pch=16, cex=1/2, bty='l')
legend("topright", c("activated", "resting"), pch=16, col=c("orange", "darkseagreen3"), bty="n")

Batch correction with Harmony

library(harmony)

## Loading required package: Rcpp

harmony_embeddings <- HarmonyMatrix(logFracs, meta_data=grpHBC, do_pca=TRUE)

## Harmony 1/10

## Harmony 2/10

## Harmony converged after 2 iterations

umHarmony <- uwot::umap(harmony_embeddings)

plot(umHarmony, col=alpha(c("orange", "darkseagreen3")[as.numeric(grpHBC)], .75), pch=16, cex=1/2, bty='l')
legend("topright", c("activated", "resting"), pch=16, col=c("orange", "darkseagreen3"), bty="n")

Integration: Processing and dimensionality reduction of each dataset separately

library(Signac)
library(Seurat)

## 
## Attaching package: 'Seurat'

## The following object is masked from 'package:SummarizedExperiment':
## 
##     Assays

library(patchwork)

## 
## Attaching package: 'patchwork'

## The following object is masked from 'package:cowplot':
## 
##     align_plots

library(ArchR)

## Loading required package: data.table

## 
## Attaching package: 'data.table'

## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

## The following object is masked from 'package:purrr':
## 
##     transpose

## The following object is masked from 'package:SummarizedExperiment':
## 
##     shift

## The following object is masked from 'package:GenomicRanges':
## 
##     shift

## The following object is masked from 'package:IRanges':
## 
##     shift

## The following objects are masked from 'package:S4Vectors':
## 
##     first, second

## Loading required package: rhdf5

## Loading required package: magrittr

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:purrr':
## 
##     set_names

## The following object is masked from 'package:tidyr':
## 
##     extract

set.seed(1234)
addArchRGenome("mm10")

## Setting default genome to Mm10.

# scATAC-seq importing and preprocessing
# from 20200610_scATAC_injury_archr_samples_pairwiseHBC_v3.Rmd
archProj <- loadArchRProject("../data/scATAC/archR_samples_pairwise_v2")

## Successfully loaded ArchRProject!

## 
##                                                    / |
##                                                  /    \
##             .                                  /      |.
##             \\\                              /        |.
##               \\\                          /           `|.
##                 \\\                      /              |.
##                   \                    /                |\
##                   \\#####\           /                  ||
##                 ==###########>      /                   ||
##                  \\##==......\    /                     ||
##             ______ =       =|__ /__                     ||      \\\
##         ,--' ,----`-,__ ___/'  --,-`-===================##========>
##        \               '        ##_______ _____ ,--,__,=##,__   ///
##         ,    __==    ___,-,__,--'#'  ==='      `-'    | ##,-/
##         -,____,---'       \\####\\________________,--\\_##,/
##            ___      .______        ______  __    __  .______      
##           /   \     |   _  \      /      ||  |  |  | |   _  \     
##          /  ^  \    |  |_)  |    |  ,----'|  |__|  | |  |_)  |    
##         /  /_\  \   |      /     |  |     |   __   | |      /     
##        /  _____  \  |  |\  \\___ |  `----.|  |  |  | |  |\  \\___.
##       /__/     \__\ | _| `._____| \______||__|  |__| | _| `._____|
## 

clHBCMerged <- archProj$clHBCMerged
plotEmbedding(ArchRProj = archProj, colorBy = "cellColData", 
              name = "clHBCMerged", embedding = "UMAP_cor")

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-4f7e696fefec-Date-2021-01-29_Time-16-32-26.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = cellColData

## Plotting Embedding

## 1 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-4f7e696fefec-Date-2021-01-29_Time-16-32-26.log

archrUmap <- ArchR::getEmbedding(archProj, embedding = "UMAP_cor")
colnames(archrUmap) <- c("UMAP1", "UMAP2")
plotEmbedding(ArchRProj = archProj, colorBy = "cellColData", 
              name = "treat", embedding = "UMAP_cor")

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-4f7e33f0a508-Date-2021-01-29_Time-16-32-29.log
## If there is an issue, please report to github with logFile!
## Getting UMAP Embedding
## ColorBy = cellColData
## Plotting Embedding
## 1 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-4f7e33f0a508-Date-2021-01-29_Time-16-32-29.log

clHBCMerged2 <- as.character(getCellColData(archProj, "clHBCMerged")$clHBCMerged)
clHBCMerged2[clHBCMerged2 == "C1_Inj"] <- "Activated"
clHBCMerged2[clHBCMerged2 == "C2_Hybrid"] <- "Hybrid"
clHBCMerged2[clHBCMerged2 == "C3_UI"] <- "Resting"
archrUmap$clHBCMerged2 <- factor(clHBCMerged2)
archrUmap$treatment <- getCellColData(archProj, "treat")$treat
peaks <- archProj@peakSet
names(peaks) <- paste0("peak",1:length(peaks))
atacPeakMat <- readRDS("../data/scATAC/peakMatArchrHbc_v2.rds")
rownames(atacPeakMat) <- names(peaks)
geneActivity <- readRDS("../data/scATAC/geneMatArchrHbc.rds")
atac <- Seurat::CreateSeuratObject(counts = atacPeakMat,
                           assay = "ATAC",
                           project = "10x_ATAC")
atac[["ACTIVITY"]] <- CreateAssayObject(counts = geneActivity)
atac$tech <- "atac"
atac$treatment <- archProj@cellColData$treat
atac$clHBCMerged <- clHBCMerged
# process gene activity
DefaultAssay(atac) <- "ACTIVITY"
atac <- FindVariableFeatures(atac)
atac <- NormalizeData(atac)
atac <- ScaleData(atac)

## Centering and scaling data matrix

# process peak matrix
DefaultAssay(atac) <- "ATAC"
# VariableFeatures(atac) <- names(which(Matrix::rowSums(atac) > 50))
# atac <- RunLSI(atac, n = 30, scale.max = NULL)
# atac <- RunUMAP(atac, reduction = "lsi", dims = 1:30)
atac <- FindTopFeatures(atac, min.cutoff = 'q0')
# Seurat 3
atac <- RunLSI(atac, n = 30) 

## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac

## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac

## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac

## Performing TF-IDF normalization
## Running SVD on TF-IDF matrix
## Scaling cell embeddings

# Seurat 4
#atac <- RunTFIDF(atac)
#atac <- RunSVD(atac, n = 30, reduction.key = "lsi_", reduction.name = "lsi")
# drop first dim as often correlated with seq depth
atac <- RunUMAP(object = atac, reduction = 'lsi', dims = 2:30)

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 16:34:01 UMAP embedding parameters a = 0.9922 b = 1.112
## 16:34:01 Read 4076 rows and found 29 numeric columns
## 16:34:01 Using Annoy for neighbor search, n_neighbors = 30
## 16:34:01 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 16:34:02 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpQk0sz5/file4f7e7be2cadb
## 16:34:02 Searching Annoy index using 1 thread, search_k = 3000
## 16:34:03 Annoy recall = 100%
## 16:34:05 Commencing smooth kNN distance calibration using 1 thread
## 16:34:07 Initializing from normalized Laplacian + noise
## 16:34:07 Commencing optimization for 500 epochs, with 152690 positive edges
## 16:34:14 Optimization finished

# scRNA-seq import
rna <- CreateSeuratObject(counts=countHBC)
rna$tech <- '10x'
rna$celltype <- grpHBC
rna <- FindVariableFeatures(rna, selection.method = "vst", nfeatures = 2000)
rna <- ScaleData(rna)

## Centering and scaling data matrix

rna <- RunPCA(rna, npcs=20)

## PC_ 1 
## Positive:  Slc6a6, Rps27, Junb, Aqp3, Btg2, Fos, Fxyd3, Apoe, Krt5, Jun 
##     Gclm, Socs3, Ier2, Gpm6a, Ubc, Ctsl, Cyr61, Ifitm3, Laptm4a, Atf3 
##     Gstm2, Csrp1, Ptn, Pltp, Fosb, Dusp1, Igfbp2, Klf4, Krt14, Gas6 
## Negative:  Tmsb10, Cystm1, Atf5, Plaur, Adam8, Krt6a, Map1b, Tuba1a, Tubb2b, Cst6 
##     Ptma, Ecm1, Stmn1, Emp3, Uchl1, Upk3bl, Ceacam1, Clcf1, Tm4sf1, Gprc5a 
##     Cldn7, Tubb6, Lypd3, Cldn3, Lamc2, Phlda2, Ncam1, Cldn4, Mal, Sult2b1 
## PC_ 2 
## Positive:  Igfbp2, Cald1, Tmsb10, Ogn, Tmsb4x, Igfbp4, Adam8, Krt5, Apoe, Plaur 
##     Fam101b, Tspan6, Abca4, Sparc, Tm4sf1, Col23a1, Lamc2, Krt6a, Clcf1, Tmem108 
##     Slc29a1, Krt14, Tubb6, Ngf, Cst6, Oaf, Tubb2b, Ecm1, Hmga2, Car12 
## Negative:  Prdx6, Adh7, Gsta3, Sftpd, Cyp2f2, Ly6a, Tst, Slc16a11, Ifitm1, Reg3g 
##     Muc1, Por, Cbr2, Adh1, Lrrc26, Selenbp1, Mgst1, Gpx2, Sult1d1, Gstp1 
##     Cbr3, Cd14, Cyp4b1, Cxcl17, Gsto1, Serpinb11, Elf3, Ly6c1, Serpina3n, Epas1 
## PC_ 3 
## Positive:  Adh7, Aldh3a1, Selenbp1, Adh1, Ly6a, Oat, Defb1, Fmo2, Cbr3, Gstp1 
##     Tmem176a, Atp1b1, Tmem176b, Efemp1, Igfbp4, Upk1b, Sparc, Ly6c1, Tgm2, Igfbp2 
##     Mgp, Kcnj16, Sult1d1, Gmpr, Cdo1, Krt15, Nrtn, Cyp4b1, Muc1, Bsg 
## Negative:  Cldn4, X1810011O10Rik, Klf4, Sfn, Zfp36, S100a10, Mast4, Atf3, Tagln2, Fosb 
##     Avpi1, Arc, Nr4a1, Cldn7, Nfat5, Egr1, Klf6, Cdkn1a, Phlda1, Lmna 
##     Crym, Neat1, Anxa1, Krt18, S100a6, Acsm4, Lypd3, Papss2, Ier2, Krt23 
## PC_ 4 
## Positive:  mt.Co1, mt.Co3, mt.Atp6, Hpgd, Pla2r1, X1500009L16Rik, Acsm4, Rps27rt, Car12, Dapl1 
##     Myo6, Fth1, Papss2, Lypd2, Insl6, Cyb5r3, Scd2, Gas6, Trp63, Nmb 
##     Eya1, Piezo2, Limd2, Crym, Sema5a, Serpinb1b, Ftl1, X4631405K08Rik, Fndc1, Scgb1c1 
## Negative:  Actg1, Anxa3, Ier3, Nfkbia, F3, Tmem176b, Krt7, Sfn, Cyr61, Cldn4 
##     Krt14, Tacstd2, Tuba1c, Tmem176a, Fos, Bhlhe40, Junb, Tnfrsf12a, Gmpr, Cish 
##     Fst, Cxcl16, Plaur, Krt17, Ubc, Trib1, Oat, Atf3, Adam8, Dusp5 
## PC_ 5 
## Positive:  Fosb, Egr1, Neat1, Irs2, Pdcd4, Nr4a1, Lmna, Klf9, Cbx3, Dst 
##     Ccdc3, Mast4, Sdc4, Timp3, Nfia, Kdm6b, Birc5, Maff, Spry2, Il6st 
##     Stmn1, Hist1h2ap, Nfat5, Sik1, Arid5a, Gm26532, Pbk, Gclc, Tpm1, Fam101b 
## Negative:  Aldh2, Gstm1, Prdx1, Aqp3, Ces1d, Tst, Gsn, Clu, Ephx1, Cpe 
##     Ifitm3, Ubc, Krt18, Prdx6, Hspb1, Galm, Fmo6, Hspa5, Ppa1, Krt5 
##     Krt8, Pltp, Mgst1, Ctsl, Cyp2a5, Alas1, Cryl1, Calml4, Gstm2, Actg1

rna <- RunUMAP(rna, dims = 1:20, min_dist=.1)

## Warning: The following arguments are not used: min_dist

## 16:34:18 UMAP embedding parameters a = 0.9922 b = 1.112
## 16:34:18 Read 3064 rows and found 20 numeric columns
## 16:34:18 Using Annoy for neighbor search, n_neighbors = 30
## 16:34:18 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 16:34:19 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpQk0sz5/file4f7e66742f67
## 16:34:19 Searching Annoy index using 1 thread, search_k = 3000
## 16:34:19 Annoy recall = 100%
## 16:34:21 Commencing smooth kNN distance calibration using 1 thread
## 16:34:22 Initializing from normalized Laplacian + noise
## 16:34:22 Commencing optimization for 500 epochs, with 122146 positive edges
## 16:34:28 Optimization finished

p1 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=TRUE) + 
  NoLegend() + 
  ggtitle("scATAC-seq")

## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.

p2 <- DimPlot(rna, group.by = "celltype", label = TRUE, repel = TRUE) + NoLegend() + ggtitle("scRNA-seq")
p1 + p2

# Markers for activated/resting HBC in scRNA-seq
FeaturePlot(rna, features = "Krtdap")

FeaturePlot(rna, features = "Krt6a")

FeaturePlot(rna, features = "Krt5")

FeaturePlot(rna, features = "Krt14")

FeaturePlot(rna, features = "Trp63")

# Markers for activated/resting HBC in scATAC-seq
FeaturePlot(atac, features = "Krtdap") + xlim(c(-5,7))

## Warning: Could not find Krtdap in the default search locations, found in
## ACTIVITY assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 14 rows containing missing values (geom_point).

FeaturePlot(atac, features = "Krt6a") + xlim(c(-5,7))

## Warning: Could not find Krt6a in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 14 rows containing missing values (geom_point).

FeaturePlot(atac, features = "Krt5") + xlim(c(-6,6))

## Warning: Could not find Krt5 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

FeaturePlot(atac, features = "Krt14") + xlim(c(-6,6))

## Warning: Could not find Krt14 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

FeaturePlot(atac, features = "Trp63") + xlim(c(-6,6))

## Warning: Could not find Trp63 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

Seurat reduced dimension plots

# RNA
pRNA1 <- DimPlot(rna, group.by = "celltype", label = FALSE, repel = TRUE) + 
  ggtitle("scRNA-seq: 24h HBC* and rHBC cells") +
  scale_color_manual(labels = c("Activated", "Regenerated\n(resting)"), 
                     values = c("steelblue", "salmon")) +
  labs(color = "Cell state")
# ggsave("../plots/scRNA24hHBC.png", width=9)

# ATAC based on Seurat
pATAC1 <- DimPlot(atac, reduction = "umap", group.by="clHBCMerged", label=FALSE) + 
  ggtitle("scATAC-seq HBC cell states") + 
  xlim(c(-6, 7))
pATAC2 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=FALSE, pt.size = .2) +
  scale_color_manual(labels = c("Injured", "Uninjured"),
                     values = c("orange", "darkseagreen3")) +
  ggtitle("scATAC-seq HBC treatment") + 
  xlim(c(-6, 7))
cowplot::plot_grid(pATAC1, pATAC2)

## Warning: Removed 14 rows containing missing values (geom_point).

## Warning: Removed 14 rows containing missing values (geom_point).

# ggsave("../plots/ATACClusters_seurat.png", width=10)

Data integration

Cell type transfer

# find anchor features
transfer.anchors <- FindTransferAnchors(reference = rna, query = atac, 
                                        features = VariableFeatures(object = rna),
                                        reference.assay = "RNA", 
                                        query.assay = "ACTIVITY", 
                                        reduction = "cca")

## Warning in RunCCA.Seurat(object1 = reference, object2 = query, features =
## features, : Running CCA on different assays

## Running CCA

## Merging objects

## Finding neighborhoods

## Finding anchors

##  Found 10751 anchors

## Filtering anchors

##  Retained 1185 anchors

## Extracting within-dataset neighbors

# predict ATAC cell type based on RNA labels
celltype.predictions <- TransferData(anchorset = transfer.anchors, 
                                     refdata = rna$celltype, 
                                     weight.reduction = atac[["lsi"]])

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

# add predicted labels to ATAC
atac <- AddMetaData(atac, metadata = celltype.predictions)
hist(atac$prediction.score.max)

atac$predicted.id <- factor(atac$predicted.id)  # to make the colors match

Data integration plots in Seurat space

pat1 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=TRUE) + 
  NoLegend() + 
  ggtitle("scATAC-seq truth") +
  xlim(c(-6,7))
pat2 <- DimPlot(atac, group.by = "predicted.id", label = TRUE, repel = TRUE) + 
  ggtitle("scATAC-seq predicted") + 
  NoLegend() + 
  scale_colour_hue(drop = FALSE) +
  xlim(c(-6,7))
pat3 <- DimPlot(atac, group.by = "clHBCMerged", label = TRUE, repel = TRUE) + 
  ggtitle("scATAC-seq ArchR clusters") + 
  NoLegend() + 
  scale_colour_hue(drop = FALSE) +
  xlim(c(-6,7))
prna4 <- DimPlot(rna, group.by = "celltype", label = TRUE, repel = TRUE) + ggtitle("scRNA-seq") + 
    NoLegend()
pat1 + pat2 + pat3 + prna4

## Warning: Removed 14 rows containing missing values (geom_point).

## Warning: Removed 14 rows containing missing values (geom_point).

## Warning: Removed 14 rows containing missing values (geom_point).

Data integration plots in ArchR space

# prediction in ArchR DR space
archrUmap$predicted.id <- celltype.predictions$predicted.id
pat1_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=treatment)) +
  geom_point(size = .75) +
  ggtitle("scATAC-seq treatment") +
  theme_classic() +
  scale_color_manual(values = c("orange", "darkseagreen3"))
pat2_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=predicted.id)) +
  geom_point(size = .75) +
  ggtitle("scATAC-seq predicted") +
  theme_classic()
pat3_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=clHBCMerged)) +
  geom_point(size = .75) +
  ggtitle("scATAC-seq clusters") +
  theme_classic()
prna1_2 <- DimPlot(rna, group.by = "celltype", label = TRUE, repel = TRUE) + 
  ggtitle("scRNA-seq") + 
  NoLegend()
pat1_2 + pat2_2 + pat3_2 + prna1_2

Joint dimensionality reduction

genes.use <- VariableFeatures(rna)
refdata <- GetAssayData(rna, assay = "RNA", slot = "data")[genes.use, ]

# refdata (input) contains a scRNA-seq expression matrix for the scRNA-seq cells.  imputation
# (output) will contain an imputed scRNA-seq matrix for each of the ATAC cells
imputation <- TransferData(anchorset = transfer.anchors, 
                           refdata = refdata, 
                           weight.reduction = atac[["lsi"]])

## Finding integration vectors

## Finding integration vector weights

## Transfering 2000 features onto reference data

# this line adds the imputed data matrix to the pbmc.atac object
atacJoint <- atac
atacJoint[["RNA"]] <- imputation
coembed <- merge(x = rna, y = atacJoint)

# Finally, we run PCA and UMAP on this combined object, to visualize the co-embedding of both
# datasets
coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE)

## Centering data matrix

coembed <- RunPCA(coembed, features = genes.use, verbose = FALSE)
coembed <- RunUMAP(coembed, dims = 1:30)

## 16:36:36 UMAP embedding parameters a = 0.9922 b = 1.112

## 16:36:36 Read 7140 rows and found 30 numeric columns

## 16:36:36 Using Annoy for neighbor search, n_neighbors = 30

## 16:36:36 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 16:36:38 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpQk0sz5/file4f7e71fba566
## 16:36:38 Searching Annoy index using 1 thread, search_k = 3000
## 16:36:42 Annoy recall = 100%
## 16:36:44 Commencing smooth kNN distance calibration using 1 thread
## 16:36:46 Initializing from normalized Laplacian + noise
## 16:36:46 Commencing optimization for 500 epochs, with 318608 positive edges
## 16:36:58 Optimization finished

coembed$celltype <- ifelse(!is.na(coembed$celltype), coembed$celltype, coembed$predicted.id)

p1 <- DimPlot(coembed, group.by = "tech")
p2 <- DimPlot(coembed, group.by = "celltype", label = TRUE, repel = TRUE)
p1 + p2

Compare predictions with truth and ArchR clusters

library(ggalluvial)

dfCellWide <- data.frame(treatment=atac$treatment,
                     predicted=atac$predicted.id,
                     cluster= archProj@cellColData$clHBCMerged)

pAlluv <- ggplot(dfCellWide,
       aes( axis1 = cluster, axis2 = predicted, axis3=treatment)) +
  geom_alluvium(aes(fill=treatment), width = 1/12) +
  scale_x_discrete(limits = c("Cluster", "Predicted", "Treatment"), expand = c(.05, .05)) +
  geom_stratum(width = 1/12, fill = "white", color = "grey") +
  geom_label(stat = "stratum", aes(label = after_stat(stratum))) +
  theme(panel.background = element_rect(fill="white"))
pAlluv

Check gene expression of chromatin-responsive genes upon injury

Here, we check the gene expression for genes that were found to be DA in the scATAC-seq data when testing for the injury effect for both the injured vs uninjured cluster and injury vs uninjured cells in the hybrid cluster.

Chromatin more accessible in injury

chromUpInjury <- read.table("../data/scATAC/sharedUpInjury.txt",
                            stringsAsFactors = FALSE)[,1]
chromUpInjury <- chromUpInjury[chromUpInjury %in% rownames(countHBC)]


plotByExpr <- function(rna, genes, counts, name="feature", scale=TRUE){
  genes <- genes[genes %in% rownames(countHBC)]
  if(scale){
    scaledCounts <- t(scale(t(counts[genes,]+1e-5)))
    scaledCounts <- scaledCounts[!is.na(rowVars(scaledCounts)),]
    yhat <- colSums(scaledCounts)
  } else {
    yhat <- colSums(counts[genes,]+1e-5)
  }
  rna <- AddMetaData(rna, yhat, col.name = name)
  FeaturePlot(rna, features = name) + theme(legend.position = "none")
}

plotByExpr(rna, chromUpInjury, countHBC, "chromUpInjurySum", scale=FALSE)

plotByExpr(rna, chromUpInjury, countHBC, "chromUpInjurySum_scaled", scale=TRUE)

Chromatin less accessible in injury

chromDownInjury <- read.table("../data/scATAC/sharedDownInjury.txt",
                            stringsAsFactors = FALSE)[,1]
chromDownInjury <- chromDownInjury[chromDownInjury %in% rownames(countHBC)]

plotByExpr(rna, chromDownInjury, countHBC, "chromDownInjurySum", scale=FALSE)

plotByExpr(rna, chromDownInjury, countHBC, "chromDownjurySum_scaled", scale=TRUE)

Check gene expression of chromatin-responsive genes in hybrid vs uninjured ATAC-seq clusters

chromHybrid <- read.table("../data/scATAC/sharedHybridGenes.txt",
                            stringsAsFactors = FALSE)[,1]
chromHybrid <- chromHybrid[chromHybrid %in% rownames(countHBC)]

plotByExpr(rna, chromHybrid, countHBC, "hybridSum", scale=FALSE)

plotByExpr(rna, chromHybrid, countHBC, "chybridSum_scaled", scale=TRUE)

ATAC-seq cluster markers

# atacMarkers <- readRDS("../../../scATAC/data/markerList_atac_ArchR_hbcClusters.rds")
atacMarkers <- readRDS("../data/scATAC/markerList_atac_ArchR_hbcClusters.rds")


pc1 <- plotByExpr(rna, atacMarkers$C1_Inj$name, countHBC, "C1_Inj", scale=TRUE)
pc2 <- plotByExpr(rna, atacMarkers$C2_Hybrid$name, countHBC, "C2_Hybrid", scale=TRUE)
pc3 <- plotByExpr(rna, atacMarkers$C3_UI$name, countHBC, "C3_UI", scale=TRUE)

pMarkersATAC <- cowplot::plot_grid(plotlist=list(pc1,pc2,pc3), nrow=1)
# ggsave("../plots/ATACMarkersOnRNA.png", width=19)

## are the ATAC markers driven by both injured and uninjured cells in hybrid cluster? Yes.
activityHybridMarkers <- as.matrix(atac[["ACTIVITY"]]@counts)[atacMarkers$C2_Hybrid$name,]
df <- data.frame(y=colSums(t(scale(t(activityHybridMarkers)))),
                 state=factor(interaction(atac$clHBCMerged, atac$treatment)))
ggplot(df, aes(x=state, y=y)) +
  geom_boxplot() +
  xlab("Group (cell state x treatment)") +
  ylab("Scaled activity of hybrid markers") +
  ggtitle("Cell state expression of hybrid markers") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, vjust=0.5))

# ggsave("../plots/atacHybridMarkerCellContribution.png", width=9)

Redo cell label prediction when assigning hybrid labels

For reference, old plot

# Seurat
pat1 + pat2 + pat3 + prna4

## Warning: Removed 14 rows containing missing values (geom_point).

## Warning: Removed 14 rows containing missing values (geom_point).

## Warning: Removed 14 rows containing missing values (geom_point).

# ArchR DR space
pat1_2 + pat2_2 + pat3_2 + prna1_2

Per hybrid cluster separately

ctHybridSub <- as.character(rna$celltype)
hlpid <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>-1)
ctHybridSub[hlpid] <- NA
hlp2id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>-1 & rna@reductions$umap@cell.embeddings[,2]<4)
ctHybridSub[hlp2id] <- "hybrid1"
hlp1id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>4)
ctHybridSub[hlp1id] <- "hybrid2"
rna$ctHybridSub <- ctHybridSub
# ctHybridSub <- as.character(rna$celltype)
# hlpid <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2] < 2)
# ctHybridSub[hlpid] <- NA
# hlp2id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2] > -1 & rna@reductions$umap@cell.embeddings[,2]<2)
# ctHybridSub[hlp2id] <- "hybrid1"
# hlp1id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2] < -4)
# ctHybridSub[hlp1id] <- "hybrid2"
# rna$ctHybridSub <- ctHybridSub


DimPlot(rna, group.by = "ctHybridSub", label = FALSE, repel = TRUE) +
  ggtitle("scRNA-seq") 

# predict ATAC cell type based on RNA labels
celltype.predictions3 <- TransferData(anchorset = transfer.anchors, 
                                     refdata = rna$ctHybridSub, 
                                     weight.reduction = atac[["lsi"]])

## Finding integration vectors

## Finding integration vector weights

## Predicting cell labels

# add predicted labels to ATAC
atac$predictedCtHybrid <- celltype.predictions3$predicted.id
atac$predScore <- celltype.predictions3$prediction.score.max
atac$predicted.id <- factor(atac$predicted.id)  # to make the colors match


# prediction in Seurat RD space
ppred1 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=FALSE) + 
  ggtitle("scATAC-seq treatment") +
  xlim(c(-6, 7))
ppred2 <- DimPlot(atac, group.by = "predictedCtHybrid", label = FALSE, repel = TRUE) + 
  ggtitle("scATAC-seq predicted state") + 
  scale_colour_hue(drop = FALSE) +
  xlim(c(-6, 7))
ppred3 <- DimPlot(atac, group.by = "clHBCMerged", label = FALSE, repel = TRUE) + 
  ggtitle("scATAC-seq cell state") + 
  scale_colour_hue(drop = FALSE) +
  xlim(c(-6, 7))
ppred4 <- DimPlot(rna, group.by = "ctHybridSub", label = FALSE, repel = TRUE) + 
  ggtitle("scRNA-seq")
ppred1 + ppred2 + ppred3 + ppred4

## Warning: Removed 14 rows containing missing values (geom_point).

## Warning: Removed 14 rows containing missing values (geom_point).

## Warning: Removed 14 rows containing missing values (geom_point).

# ggsave("../plots/cellTransferSubHybrid.png", width=10)

## size of cell is prediction prob
DimPlot(atac, group.by = "predictedCtHybrid", label = TRUE, repel = TRUE, pt.size=celltype.predictions3$prediction.score.max*1.25) + ggtitle("scATAC-seq predicted") + 
    NoLegend() + scale_colour_hue(drop = FALSE) + xlim(c(-6,6))

## Warning: Removed 55 rows containing missing values (geom_point).

boxplot(celltype.predictions3$prediction.score.max ~ celltype.predictions3$predicted.id,
        xlab="Predicted state", ylab="Probability")

# prediction in ArchR DR space
archrUmap$predicted.id <- celltype.predictions3$predicted.id
ppred1_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=treatment)) +
  geom_point(size = .75) +
  ggtitle("scATAC-seq treatment") +
  theme_classic() +
  scale_color_manual(values = c("orange", "darkseagreen3"))
ppred2_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=predicted.id)) +
  geom_point(size = .5) +
  ggtitle("scATAC-seq predicted") +
  theme_classic()
ppred3_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=clHBCMerged)) +
  geom_point(size = .75) +
  ggtitle("scATAC-seq clusters") +
  theme_classic()
ppred4_2 <- DimPlot(rna, group.by = "ctHybridSub", label = FALSE, repel = TRUE, pt.size=.1) +
  ggtitle("scRNA-seq") +
  theme(axis.title = element_text(size = 11),
        axis.text = element_text(size = 9),
        plot.title = element_text(face = "plain", size=13))
ppred1_2 + ppred2_2 + ppred3_2 + ppred4_2

ggsave("../plots/predictionsRNAATAC_archRSpace.png", width=12, height=9)

## size of cell is prediction prob
 ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=predicted.id)) +
  geom_point(size = celltype.predictions3$prediction.score.max*1.25) +
  ggtitle("scATAC-seq predicted") +
  theme_classic()

boxplot(celltype.predictions3$prediction.score.max ~ celltype.predictions3$predicted.id,
        xlab="Predicted state", ylab="Probability")

## distribution of injured and uninjured cells in each cluster
df <- data.frame(treatment=atac$treatment,
                 cluster=atac$predictedCtHybrid)
dfTotal <- df %>% 
  group_by(cluster) %>%
  mutate(totalCells=n()) %>%
  group_by(cluster, treatment) %>%
  summarize(fracCells = n() / unique(totalCells))

## `summarise()` has grouped output by 'cluster'. You can override using the `.groups` argument.

ggplot(dfTotal, aes(x=treatment, fill=treatment, y=fracCells)) + 
  geom_bar(position="dodge", stat="identity") +
  scale_fill_manual(values = c("orange", "darkseagreen3")) +
  facet_wrap(.~cluster) +
  ylab("Fraction of cells") +
  theme_classic()

Joint dimensionality reduction

atacJoint <- atac
DefaultAssay(atacJoint) <- "ATAC"
VariableFeatures(atacJoint) <- names(which(Matrix::rowSums(atacJoint) > 100))
atacJoint <- RunLSI(atacJoint, n = 50, scale.max = NULL)

## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac

## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac

## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac

## Performing TF-IDF normalization

## Running SVD on TF-IDF matrix

## Scaling cell embeddings

atacJoint <- RunUMAP(atacJoint, reduction = "lsi", dims = 2:50)

## 16:41:20 UMAP embedding parameters a = 0.9922 b = 1.112

## 16:41:20 Read 4076 rows and found 49 numeric columns

## 16:41:20 Using Annoy for neighbor search, n_neighbors = 30

## 16:41:21 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 16:41:22 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpQk0sz5/file4f7e1979a4a5
## 16:41:22 Searching Annoy index using 1 thread, search_k = 3000
## 16:41:23 Annoy recall = 100%
## 16:41:28 Commencing smooth kNN distance calibration using 1 thread
## 16:41:30 Initializing from normalized Laplacian + noise
## 16:41:30 Commencing optimization for 500 epochs, with 155968 positive edges
## 16:41:36 Optimization finished

genes.use <- VariableFeatures(rna)
refdata <- GetAssayData(rna, assay = "RNA", slot = "data")[genes.use, ]

# refdata (input) contains a scRNA-seq expression matrix for the scRNA-seq cells.  imputation
# (output) will contain an imputed scRNA-seq matrix for each of the ATAC cells
imputation <- TransferData(anchorset = transfer.anchors, 
                           refdata = refdata, 
                           weight.reduction = atacJoint[["lsi"]])

## Finding integration vectors
## Finding integration vector weights
## Transfering 2000 features onto reference data

# this line adds the imputed data matrix to the pbmc.atac object
atacJoint[["RNA"]] <- imputation
coembed <- merge(x = rna, y = atacJoint)

# Finally, we run PCA and UMAP on this combined object, to visualize the co-embedding of both
# datasets
coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE)

## Centering data matrix

coembed <- RunPCA(coembed, features = genes.use, verbose = FALSE)
coembed <- RunUMAP(coembed, dims = 1:30)

## 16:42:16 UMAP embedding parameters a = 0.9922 b = 1.112
## 16:42:16 Read 7140 rows and found 30 numeric columns
## 16:42:16 Using Annoy for neighbor search, n_neighbors = 30
## 16:42:16 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 16:42:18 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpQk0sz5/file4f7e12859b4
## 16:42:18 Searching Annoy index using 1 thread, search_k = 3000
## 16:42:20 Annoy recall = 100%
## 16:42:23 Commencing smooth kNN distance calibration using 1 thread
## 16:42:25 Initializing from normalized Laplacian + noise
## 16:42:25 Commencing optimization for 500 epochs, with 323034 positive edges
## 16:42:37 Optimization finished

celltype.predictionsJoint <- TransferData(anchorset = transfer.anchors, 
                                     refdata = rna$ctHybridSub, 
                                     weight.reduction = atacJoint[["lsi"]],
                                     dims=1:50)

## Finding integration vectors
## Finding integration vector weights
## Predicting cell labels

# fill in predicted labels: using joint DR or procedure above is very similar.
coembed$ctHybridSubPred <- coembed$ctHybridSub
## based on joint DR
# coembed@meta.data[rownames(celltype.predictionsJoint), "ctHybridSubPred"] <- celltype.predictionsJoint$predicted.id
## ones we had before
coembed@meta.data[rownames(celltype.predictions3), "ctHybridSubPred"] <- celltype.predictions3$predicted.id


p1 <- DimPlot(coembed, group.by = "tech") + ggtitle("Modality")
p2 <- DimPlot(coembed, group.by = "ctHybridSub", label = TRUE, repel = TRUE) + ggtitle("RNA clusters")
p3 <- DimPlot(coembed, group.by = "clHBCMerged", label = TRUE, repel = TRUE) + ggtitle("ATAC clusters")
p4 <- DimPlot(coembed, group.by = "ctHybridSubPred", label = TRUE, repel = TRUE) + ggtitle("Transferred clusters")
p1 + p2 + p3 + p4

ggsave("../plots/jointDimRed_v1.png", width=10, height=6)

Markers for sub-hybrid clusters

## add annotation according to prediction
predictedAnn <- archProj@cellColData$clHBCMerged
predictedAnn[predictedAnn == "C3_UI"] <- "resting"
predictedAnn[predictedAnn == "C1_inj"] <- "activated"
hybId <- predictedAnn == "C2_hybrid"
predictedAnn[which(hybId)] <- NA
predictedAnn[hybId & 
               !celltype.predictions3$predicted.id == "hybrid1" &
               !celltype.predictions3$predicted.id == "hybrid2"] <- "hybrid"
predictedAnn[which(celltype.predictions3$predicted.id == "hybrid1")] <- "hybrid1"
predictedAnn[which(celltype.predictions3$predicted.id == "hybrid2")] <- "hybrid2"
table(predictedAnn)

## predictedAnn
##    C1_Inj C2_Hybrid   hybrid1   hybrid2   resting 
##      1485       933       197       223      1238

## visualize with Seurat
atac$predAnn <- predictedAnn
DimPlot(atac, group.by = "predAnn", label = TRUE, repel = TRUE) + ggtitle("scATAC-seq predicted") + 
    NoLegend() + scale_colour_hue(drop = FALSE) + xlim(c(-6,7))

## Warning: Removed 14 rows containing missing values (geom_point).

## distribution of injured and uninjured cells in each cluster
df <- data.frame(treatment=atac$treatment,
                 cluster=atac$predAnn)
dfTotal <- df %>% 
  group_by(cluster) %>%
  mutate(totalCells=n()) %>%
  group_by(cluster, treatment) %>%
  summarize(fracCells = n() / unique(totalCells))

## `summarise()` has grouped output by 'cluster'. You can override using the `.groups` argument.

ggplot(dfTotal, aes(x=treatment, fill=treatment, y=fracCells)) + 
    geom_bar(position="dodge", stat="identity") +
  facet_wrap(.~cluster) +
  ylab("Fraction of cells")

# saveRDS(celltype.predictions3, file="../data/scATAC/20201105_celltypePredictions3.rds")


## are the ATAC markers driven by both injured and uninjured cells in all hybrid cluster? Yes.
activityHybridMarkers <- as.matrix(atac[["ACTIVITY"]]@counts)[atacMarkers$C2_Hybrid$name,]
df <- data.frame(y=colSums(t(scale(t(activityHybridMarkers)))),
                 state=factor(interaction(atac$predAnn, atac$treatment)))
ggplot(df, aes(x=state, y=y)) +
  geom_boxplot() +
  xlab("Group (cell state x treatment)") +
  ylab("Scaled activity of hybrid markers") +
  ggtitle("Cell state expression of hybrid markers") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=45, vjust=0.5))

# ggsave("../plots/atacHybridMarkerCellContributionAllHybridClusters.png", width=9)

New peak calling based on sub-clustering of hybrid cells as suggested by RNA

archProj@cellColData$predAnn <- predictedAnn

# pseudobulk for peak calling
archProj <- addGroupCoverages(ArchRProj = archProj, groupBy = "predAnn", 
                          minReplicates = 2, maxReplicates = 6,
                          force = TRUE, 
                          verbose = FALSE)

## ArchR logging to : ArchRLogs/ArchR-addGroupCoverages-4f7e3804e61b-Date-2021-01-29_Time-16-42-46.log
## If there is an issue, please report to github with logFile!

## C1_Inj (1 of 5) : CellGroups N = 4

## C2_Hybrid (2 of 5) : CellGroups N = 6

## hybrid1 (3 of 5) : CellGroups N = 3

## hybrid2 (4 of 5) : CellGroups N = 2

## resting (5 of 5) : CellGroups N = 5

## 2021-01-29 16:42:48 : Creating Coverage Files!, 0.031 mins elapsed.

## 2021-01-29 16:42:48 : Batch Execution w/ safelapply!, 0.031 mins elapsed.

## Number of Cells = 500

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 263

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 163

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 54

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 269

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 252

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 140

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 113

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 83

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 76

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 94

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 57

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 46

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 170

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 53

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 500

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 355

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 178

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 54

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## Number of Cells = 51

## Coverage File Exists!

## Added Coverage Group

## Added Metadata Group

## Added ArrowCoverage Class

## Added Coverage/Info

## Added Coverage/Info/CellNames

## 2021-01-29 16:50:31 : Adding Kmer Bias to Coverage Files!, 7.75 mins elapsed.

## Kmer Bias chr1 (1 of 21)

## chr1 Coverage File chr1 (1 of 20)
## Coverage File chr1 (2 of 20)
## Coverage File chr1 (3 of 20)
## Coverage File chr1 (4 of 20)
## Coverage File chr1 (5 of 20)
## Coverage File chr1 (6 of 20)
## Coverage File chr1 (7 of 20)
## Coverage File chr1 (8 of 20)
## Coverage File chr1 (9 of 20)
## Coverage File chr1 (10 of 20)
## Coverage File chr1 (11 of 20)
## Coverage File chr1 (12 of 20)
## Coverage File chr1 (13 of 20)
## Coverage File chr1 (14 of 20)
## Coverage File chr1 (15 of 20)
## Coverage File chr1 (16 of 20)
## Coverage File chr1 (17 of 20)
## Coverage File chr1 (18 of 20)
## Coverage File chr1 (19 of 20)
## Coverage File chr1 (20 of 20)
## Kmer Bias chr10 (2 of 21)
## chr10 Coverage File chr10 (1 of 20)
## Coverage File chr10 (2 of 20)
## Coverage File chr10 (3 of 20)
## Coverage File chr10 (4 of 20)
## Coverage File chr10 (5 of 20)
## Coverage File chr10 (6 of 20)
## Coverage File chr10 (7 of 20)
## Coverage File chr10 (8 of 20)
## Coverage File chr10 (9 of 20)
## Coverage File chr10 (10 of 20)
## Coverage File chr10 (11 of 20)
## Coverage File chr10 (12 of 20)
## Coverage File chr10 (13 of 20)
## Coverage File chr10 (14 of 20)
## Coverage File chr10 (15 of 20)
## Coverage File chr10 (16 of 20)
## Coverage File chr10 (17 of 20)
## Coverage File chr10 (18 of 20)
## Coverage File chr10 (19 of 20)
## Coverage File chr10 (20 of 20)
## Kmer Bias chr11 (3 of 21)
## chr11 Coverage File chr11 (1 of 20)
## Coverage File chr11 (2 of 20)
## Coverage File chr11 (3 of 20)
## Coverage File chr11 (4 of 20)
## Coverage File chr11 (5 of 20)
## Coverage File chr11 (6 of 20)
## Coverage File chr11 (7 of 20)
## Coverage File chr11 (8 of 20)
## Coverage File chr11 (9 of 20)
## Coverage File chr11 (10 of 20)
## Coverage File chr11 (11 of 20)
## Coverage File chr11 (12 of 20)
## Coverage File chr11 (13 of 20)
## Coverage File chr11 (14 of 20)
## Coverage File chr11 (15 of 20)
## Coverage File chr11 (16 of 20)
## Coverage File chr11 (17 of 20)
## Coverage File chr11 (18 of 20)
## Coverage File chr11 (19 of 20)
## Coverage File chr11 (20 of 20)
## Kmer Bias chr12 (4 of 21)
## chr12 Coverage File chr12 (1 of 20)
## Coverage File chr12 (2 of 20)
## Coverage File chr12 (3 of 20)
## Coverage File chr12 (4 of 20)
## Coverage File chr12 (5 of 20)
## Coverage File chr12 (6 of 20)
## Coverage File chr12 (7 of 20)
## Coverage File chr12 (8 of 20)
## Coverage File chr12 (9 of 20)
## Coverage File chr12 (10 of 20)
## Coverage File chr12 (11 of 20)
## Coverage File chr12 (12 of 20)
## Coverage File chr12 (13 of 20)
## Coverage File chr12 (14 of 20)
## Coverage File chr12 (15 of 20)
## Coverage File chr12 (16 of 20)
## Coverage File chr12 (17 of 20)
## Coverage File chr12 (18 of 20)
## Coverage File chr12 (19 of 20)
## Coverage File chr12 (20 of 20)
## Kmer Bias chr13 (5 of 21)
## chr13 Coverage File chr13 (1 of 20)
## Coverage File chr13 (2 of 20)
## Coverage File chr13 (3 of 20)
## Coverage File chr13 (4 of 20)
## Coverage File chr13 (5 of 20)
## Coverage File chr13 (6 of 20)
## Coverage File chr13 (7 of 20)
## Coverage File chr13 (8 of 20)
## Coverage File chr13 (9 of 20)
## Coverage File chr13 (10 of 20)
## Coverage File chr13 (11 of 20)
## Coverage File chr13 (12 of 20)
## Coverage File chr13 (13 of 20)
## Coverage File chr13 (14 of 20)
## Coverage File chr13 (15 of 20)
## Coverage File chr13 (16 of 20)
## Coverage File chr13 (17 of 20)
## Coverage File chr13 (18 of 20)
## Coverage File chr13 (19 of 20)
## Coverage File chr13 (20 of 20)
## Kmer Bias chr14 (6 of 21)
## chr14 Coverage File chr14 (1 of 20)
## Coverage File chr14 (2 of 20)
## Coverage File chr14 (3 of 20)
## Coverage File chr14 (4 of 20)
## Coverage File chr14 (5 of 20)
## Coverage File chr14 (6 of 20)
## Coverage File chr14 (7 of 20)
## Coverage File chr14 (8 of 20)
## Coverage File chr14 (9 of 20)
## Coverage File chr14 (10 of 20)
## Coverage File chr14 (11 of 20)
## Coverage File chr14 (12 of 20)
## Coverage File chr14 (13 of 20)
## Coverage File chr14 (14 of 20)
## Coverage File chr14 (15 of 20)
## Coverage File chr14 (16 of 20)
## Coverage File chr14 (17 of 20)
## Coverage File chr14 (18 of 20)
## Coverage File chr14 (19 of 20)
## Coverage File chr14 (20 of 20)
## Kmer Bias chr15 (7 of 21)
## chr15 Coverage File chr15 (1 of 20)
## Coverage File chr15 (2 of 20)
## Coverage File chr15 (3 of 20)
## Coverage File chr15 (4 of 20)
## Coverage File chr15 (5 of 20)
## Coverage File chr15 (6 of 20)
## Coverage File chr15 (7 of 20)
## Coverage File chr15 (8 of 20)
## Coverage File chr15 (9 of 20)
## Coverage File chr15 (10 of 20)
## Coverage File chr15 (11 of 20)
## Coverage File chr15 (12 of 20)
## Coverage File chr15 (13 of 20)
## Coverage File chr15 (14 of 20)
## Coverage File chr15 (15 of 20)
## Coverage File chr15 (16 of 20)
## Coverage File chr15 (17 of 20)
## Coverage File chr15 (18 of 20)
## Coverage File chr15 (19 of 20)
## Coverage File chr15 (20 of 20)
## Kmer Bias chr16 (8 of 21)
## chr16 Coverage File chr16 (1 of 20)
## Coverage File chr16 (2 of 20)
## Coverage File chr16 (3 of 20)
## Coverage File chr16 (4 of 20)
## Coverage File chr16 (5 of 20)
## Coverage File chr16 (6 of 20)
## Coverage File chr16 (7 of 20)
## Coverage File chr16 (8 of 20)
## Coverage File chr16 (9 of 20)
## Coverage File chr16 (10 of 20)
## Coverage File chr16 (11 of 20)
## Coverage File chr16 (12 of 20)
## Coverage File chr16 (13 of 20)
## Coverage File chr16 (14 of 20)
## Coverage File chr16 (15 of 20)
## Coverage File chr16 (16 of 20)
## Coverage File chr16 (17 of 20)
## Coverage File chr16 (18 of 20)
## Coverage File chr16 (19 of 20)
## Coverage File chr16 (20 of 20)
## Kmer Bias chr17 (9 of 21)
## chr17 Coverage File chr17 (1 of 20)
## Coverage File chr17 (2 of 20)
## Coverage File chr17 (3 of 20)
## Coverage File chr17 (4 of 20)
## Coverage File chr17 (5 of 20)
## Coverage File chr17 (6 of 20)
## Coverage File chr17 (7 of 20)
## Coverage File chr17 (8 of 20)
## Coverage File chr17 (9 of 20)
## Coverage File chr17 (10 of 20)
## Coverage File chr17 (11 of 20)
## Coverage File chr17 (12 of 20)
## Coverage File chr17 (13 of 20)
## Coverage File chr17 (14 of 20)
## Coverage File chr17 (15 of 20)
## Coverage File chr17 (16 of 20)
## Coverage File chr17 (17 of 20)
## Coverage File chr17 (18 of 20)
## Coverage File chr17 (19 of 20)
## Coverage File chr17 (20 of 20)
## Kmer Bias chr18 (10 of 21)
## chr18 Coverage File chr18 (1 of 20)
## Coverage File chr18 (2 of 20)
## Coverage File chr18 (3 of 20)
## Coverage File chr18 (4 of 20)
## Coverage File chr18 (5 of 20)
## Coverage File chr18 (6 of 20)
## Coverage File chr18 (7 of 20)
## Coverage File chr18 (8 of 20)
## Coverage File chr18 (9 of 20)
## Coverage File chr18 (10 of 20)
## Coverage File chr18 (11 of 20)
## Coverage File chr18 (12 of 20)
## Coverage File chr18 (13 of 20)
## Coverage File chr18 (14 of 20)
## Coverage File chr18 (15 of 20)
## Coverage File chr18 (16 of 20)
## Coverage File chr18 (17 of 20)
## Coverage File chr18 (18 of 20)
## Coverage File chr18 (19 of 20)
## Coverage File chr18 (20 of 20)
## Kmer Bias chr19 (11 of 21)
## chr19 Coverage File chr19 (1 of 20)
## Coverage File chr19 (2 of 20)
## Coverage File chr19 (3 of 20)
## Coverage File chr19 (4 of 20)
## Coverage File chr19 (5 of 20)
## Coverage File chr19 (6 of 20)
## Coverage File chr19 (7 of 20)
## Coverage File chr19 (8 of 20)
## Coverage File chr19 (9 of 20)
## Coverage File chr19 (10 of 20)
## Coverage File chr19 (11 of 20)
## Coverage File chr19 (12 of 20)
## Coverage File chr19 (13 of 20)
## Coverage File chr19 (14 of 20)
## Coverage File chr19 (15 of 20)
## Coverage File chr19 (16 of 20)
## Coverage File chr19 (17 of 20)
## Coverage File chr19 (18 of 20)
## Coverage File chr19 (19 of 20)
## Coverage File chr19 (20 of 20)
## Kmer Bias chr2 (12 of 21)
## chr2 Coverage File chr2 (1 of 20)
## Coverage File chr2 (2 of 20)
## Coverage File chr2 (3 of 20)
## Coverage File chr2 (4 of 20)
## Coverage File chr2 (5 of 20)
## Coverage File chr2 (6 of 20)
## Coverage File chr2 (7 of 20)
## Coverage File chr2 (8 of 20)
## Coverage File chr2 (9 of 20)
## Coverage File chr2 (10 of 20)
## Coverage File chr2 (11 of 20)
## Coverage File chr2 (12 of 20)
## Coverage File chr2 (13 of 20)
## Coverage File chr2 (14 of 20)
## Coverage File chr2 (15 of 20)
## Coverage File chr2 (16 of 20)
## Coverage File chr2 (17 of 20)
## Coverage File chr2 (18 of 20)
## Coverage File chr2 (19 of 20)
## Coverage File chr2 (20 of 20)
## Kmer Bias chr3 (13 of 21)
## chr3 Coverage File chr3 (1 of 20)
## Coverage File chr3 (2 of 20)
## Coverage File chr3 (3 of 20)
## Coverage File chr3 (4 of 20)
## Coverage File chr3 (5 of 20)
## Coverage File chr3 (6 of 20)
## Coverage File chr3 (7 of 20)
## Coverage File chr3 (8 of 20)
## Coverage File chr3 (9 of 20)
## Coverage File chr3 (10 of 20)
## Coverage File chr3 (11 of 20)
## Coverage File chr3 (12 of 20)
## Coverage File chr3 (13 of 20)
## Coverage File chr3 (14 of 20)
## Coverage File chr3 (15 of 20)
## Coverage File chr3 (16 of 20)
## Coverage File chr3 (17 of 20)
## Coverage File chr3 (18 of 20)
## Coverage File chr3 (19 of 20)
## Coverage File chr3 (20 of 20)
## Kmer Bias chr4 (14 of 21)
## chr4 Coverage File chr4 (1 of 20)
## Coverage File chr4 (2 of 20)
## Coverage File chr4 (3 of 20)
## Coverage File chr4 (4 of 20)
## Coverage File chr4 (5 of 20)
## Coverage File chr4 (6 of 20)
## Coverage File chr4 (7 of 20)
## Coverage File chr4 (8 of 20)
## Coverage File chr4 (9 of 20)
## Coverage File chr4 (10 of 20)
## Coverage File chr4 (11 of 20)
## Coverage File chr4 (12 of 20)
## Coverage File chr4 (13 of 20)
## Coverage File chr4 (14 of 20)
## Coverage File chr4 (15 of 20)
## Coverage File chr4 (16 of 20)
## Coverage File chr4 (17 of 20)
## Coverage File chr4 (18 of 20)
## Coverage File chr4 (19 of 20)
## Coverage File chr4 (20 of 20)
## Kmer Bias chr5 (15 of 21)
## chr5 Coverage File chr5 (1 of 20)
## Coverage File chr5 (2 of 20)
## Coverage File chr5 (3 of 20)
## Coverage File chr5 (4 of 20)
## Coverage File chr5 (5 of 20)
## Coverage File chr5 (6 of 20)
## Coverage File chr5 (7 of 20)
## Coverage File chr5 (8 of 20)
## Coverage File chr5 (9 of 20)
## Coverage File chr5 (10 of 20)
## Coverage File chr5 (11 of 20)
## Coverage File chr5 (12 of 20)
## Coverage File chr5 (13 of 20)
## Coverage File chr5 (14 of 20)
## Coverage File chr5 (15 of 20)
## Coverage File chr5 (16 of 20)
## Coverage File chr5 (17 of 20)
## Coverage File chr5 (18 of 20)
## Coverage File chr5 (19 of 20)
## Coverage File chr5 (20 of 20)
## Kmer Bias chr6 (16 of 21)
## chr6 Coverage File chr6 (1 of 20)
## Coverage File chr6 (2 of 20)
## Coverage File chr6 (3 of 20)
## Coverage File chr6 (4 of 20)
## Coverage File chr6 (5 of 20)
## Coverage File chr6 (6 of 20)
## Coverage File chr6 (7 of 20)
## Coverage File chr6 (8 of 20)
## Coverage File chr6 (9 of 20)
## Coverage File chr6 (10 of 20)
## Coverage File chr6 (11 of 20)
## Coverage File chr6 (12 of 20)
## Coverage File chr6 (13 of 20)
## Coverage File chr6 (14 of 20)
## Coverage File chr6 (15 of 20)
## Coverage File chr6 (16 of 20)
## Coverage File chr6 (17 of 20)
## Coverage File chr6 (18 of 20)
## Coverage File chr6 (19 of 20)
## Coverage File chr6 (20 of 20)
## Kmer Bias chr7 (17 of 21)
## chr7 Coverage File chr7 (1 of 20)
## Coverage File chr7 (2 of 20)
## Coverage File chr7 (3 of 20)
## Coverage File chr7 (4 of 20)
## Coverage File chr7 (5 of 20)
## Coverage File chr7 (6 of 20)
## Coverage File chr7 (7 of 20)
## Coverage File chr7 (8 of 20)
## Coverage File chr7 (9 of 20)
## Coverage File chr7 (10 of 20)
## Coverage File chr7 (11 of 20)
## Coverage File chr7 (12 of 20)
## Coverage File chr7 (13 of 20)
## Coverage File chr7 (14 of 20)
## Coverage File chr7 (15 of 20)
## Coverage File chr7 (16 of 20)
## Coverage File chr7 (17 of 20)
## Coverage File chr7 (18 of 20)
## Coverage File chr7 (19 of 20)
## Coverage File chr7 (20 of 20)
## Kmer Bias chr8 (18 of 21)
## chr8 Coverage File chr8 (1 of 20)
## Coverage File chr8 (2 of 20)
## Coverage File chr8 (3 of 20)
## Coverage File chr8 (4 of 20)
## Coverage File chr8 (5 of 20)
## Coverage File chr8 (6 of 20)
## Coverage File chr8 (7 of 20)
## Coverage File chr8 (8 of 20)
## Coverage File chr8 (9 of 20)
## Coverage File chr8 (10 of 20)
## Coverage File chr8 (11 of 20)
## Coverage File chr8 (12 of 20)
## Coverage File chr8 (13 of 20)
## Coverage File chr8 (14 of 20)
## Coverage File chr8 (15 of 20)
## Coverage File chr8 (16 of 20)
## Coverage File chr8 (17 of 20)
## Coverage File chr8 (18 of 20)
## Coverage File chr8 (19 of 20)
## Coverage File chr8 (20 of 20)
## Kmer Bias chr9 (19 of 21)
## chr9 Coverage File chr9 (1 of 20)
## Coverage File chr9 (2 of 20)
## Coverage File chr9 (3 of 20)
## Coverage File chr9 (4 of 20)
## Coverage File chr9 (5 of 20)
## Coverage File chr9 (6 of 20)
## Coverage File chr9 (7 of 20)
## Coverage File chr9 (8 of 20)
## Coverage File chr9 (9 of 20)
## Coverage File chr9 (10 of 20)
## Coverage File chr9 (11 of 20)
## Coverage File chr9 (12 of 20)
## Coverage File chr9 (13 of 20)
## Coverage File chr9 (14 of 20)
## Coverage File chr9 (15 of 20)
## Coverage File chr9 (16 of 20)
## Coverage File chr9 (17 of 20)
## Coverage File chr9 (18 of 20)
## Coverage File chr9 (19 of 20)
## Coverage File chr9 (20 of 20)
## Kmer Bias chrX (20 of 21)
## chrX Coverage File chrX (1 of 20)
## Coverage File chrX (2 of 20)
## Coverage File chrX (3 of 20)
## Coverage File chrX (4 of 20)
## Coverage File chrX (5 of 20)
## Coverage File chrX (6 of 20)
## Coverage File chrX (7 of 20)
## Coverage File chrX (8 of 20)
## Coverage File chrX (9 of 20)
## Coverage File chrX (10 of 20)
## Coverage File chrX (11 of 20)
## Coverage File chrX (12 of 20)
## Coverage File chrX (13 of 20)
## Coverage File chrX (14 of 20)
## Coverage File chrX (15 of 20)
## Coverage File chrX (16 of 20)
## Coverage File chrX (17 of 20)
## Coverage File chrX (18 of 20)
## Coverage File chrX (19 of 20)
## Coverage File chrX (20 of 20)
## Kmer Bias chrY (21 of 21)
## chrY Coverage File chrY (1 of 20)
## Coverage File chrY (2 of 20)
## Coverage File chrY (3 of 20)
## Coverage File chrY (4 of 20)
## Coverage File chrY (5 of 20)
## Coverage File chrY (6 of 20)
## Coverage File chrY (7 of 20)
## Coverage File chrY (8 of 20)
## Coverage File chrY (9 of 20)
## Coverage File chrY (10 of 20)
## Coverage File chrY (11 of 20)
## Coverage File chrY (12 of 20)
## Coverage File chrY (13 of 20)
## Coverage File chrY (14 of 20)
## Coverage File chrY (15 of 20)
## Coverage File chrY (16 of 20)
## Coverage File chrY (17 of 20)
## Coverage File chrY (18 of 20)
## Coverage File chrY (19 of 20)
## Coverage File chrY (20 of 20)
## Completed Kmer Bias Calculation
## Adding Kmer Bias (1 of 20)
## Adding Kmer Bias (2 of 20)
## Adding Kmer Bias (3 of 20)
## Adding Kmer Bias (4 of 20)
## Adding Kmer Bias (5 of 20)
## Adding Kmer Bias (6 of 20)
## Adding Kmer Bias (7 of 20)
## Adding Kmer Bias (8 of 20)
## Adding Kmer Bias (9 of 20)
## Adding Kmer Bias (10 of 20)
## Adding Kmer Bias (11 of 20)
## Adding Kmer Bias (12 of 20)
## Adding Kmer Bias (13 of 20)
## Adding Kmer Bias (14 of 20)
## Adding Kmer Bias (15 of 20)
## Adding Kmer Bias (16 of 20)
## Adding Kmer Bias (17 of 20)
## Adding Kmer Bias (18 of 20)
## Adding Kmer Bias (19 of 20)
## Adding Kmer Bias (20 of 20)
## 2021-01-29 16:58:01 : Finished Creation of Coverage Files!, 15.252 mins elapsed.
## ArchR logging successful to : ArchRLogs/ArchR-addGroupCoverages-4f7e3804e61b-Date-2021-01-29_Time-16-42-46.log

# peak calling
archProj <- addReproduciblePeakSet(
    ArchRProj = archProj, 
    groupBy = "predAnn", 
    pathToMacs2 = "/Users/koenvandenberge/opt/anaconda3/bin/macs2",
    verbose = FALSE)

## ArchR logging to : ArchRLogs/ArchR-addReproduciblePeakSet-4f7e499f1358-Date-2021-01-29_Time-16-58-01.log
## If there is an issue, please report to github with logFile!
## Calling Peaks with Macs2
## 2021-01-29 16:58:02 : Batch Execution w/ safelapply!, 0 mins elapsed.
## 2021-01-29 16:58:02 : Group 1 of 20, Calling Peaks with MACS2!, 0 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.s1Inj-1 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.s1Inj-1.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 16:59:50 : Group 2 of 20, Calling Peaks with MACS2!, 1.8 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.s2Inj-2 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.s2Inj-2.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:01:26 : Group 3 of 20, Calling Peaks with MACS2!, 3.398 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.s3Inj-3 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.s3Inj-3.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:02:43 : Group 4 of 20, Calling Peaks with MACS2!, 4.691 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.Other-4 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.Other-4.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:03:34 : Group 5 of 20, Calling Peaks with MACS2!, 5.538 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s1Inj-5 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s1Inj-5.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:04:41 : Group 6 of 20, Calling Peaks with MACS2!, 6.651 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s3Inj-6 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s3Inj-6.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:06:02 : Group 7 of 20, Calling Peaks with MACS2!, 7.999 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s2UI-7 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s2UI-7.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:06:47 : Group 8 of 20, Calling Peaks with MACS2!, 8.761 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s3UI-8 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s3UI-8.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:07:32 : Group 9 of 20, Calling Peaks with MACS2!, 9.501 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s2Inj-9 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s2Inj-9.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:08:22 : Group 10 of 20, Calling Peaks with MACS2!, 10.335 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s1UI-10 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s1UI-10.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:09:15 : Group 11 of 20, Calling Peaks with MACS2!, 11.213 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid1._.s3UI-11 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid1._.s3UI-11.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:10:20 : Group 12 of 20, Calling Peaks with MACS2!, 12.304 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid1._.s3Inj-12 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid1._.s3Inj-12.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:11:20 : Group 13 of 20, Calling Peaks with MACS2!, 13.303 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid1._.Other-13 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid1._.Other-13.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:12:17 : Group 14 of 20, Calling Peaks with MACS2!, 14.255 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid2._.s3UI-14 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid2._.s3UI-14.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:13:39 : Group 15 of 20, Calling Peaks with MACS2!, 15.625 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid2._.1-15 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid2._.1-15.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:14:14 : Group 16 of 20, Calling Peaks with MACS2!, 16.21 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s2UI-16 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s2UI-16.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:15:33 : Group 17 of 20, Calling Peaks with MACS2!, 17.512 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s1UI-17 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s1UI-17.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:16:42 : Group 18 of 20, Calling Peaks with MACS2!, 18.678 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s3UI-18 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s3UI-18.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:17:44 : Group 19 of 20, Calling Peaks with MACS2!, 19.699 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s3Inj-19 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s3Inj-19.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2021-01-29 17:18:26 : Group 20 of 20, Calling Peaks with MACS2!, 20.407 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.Other-20 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.Other-20.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC

## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/C1_Inj-reproduciblePeaks.gr.rds"

## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC

## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/C2_Hybrid-reproduciblePeaks.gr.rds"

## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC

## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/hybrid1-reproduciblePeaks.gr.rds"

## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC

## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/hybrid2-reproduciblePeaks.gr.rds"

## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC

## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/resting-reproduciblePeaks.gr.rds"

## Converged after 8 iterations!
## Plotting Ggplot!

archProj <- addPeakMatrix(archProj)

## ArchR logging to : ArchRLogs/ArchR-addPeakMatrix-4f7ee62401e-Date-2021-01-29_Time-17-19-59.log
## If there is an issue, please report to github with logFile!
## 2021-01-29 17:19:59 : Batch Execution w/ safelapply!, 0 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-01-29 17:20:08 : Adding s1Inj to PeakMatrix for Chr (1 of 20)!, 0.011 mins elapsed.
## 2021-01-29 17:20:11 : Adding s1Inj to PeakMatrix for Chr (2 of 20)!, 0.057 mins elapsed.
## 2021-01-29 17:20:13 : Adding s1Inj to PeakMatrix for Chr (3 of 20)!, 0.099 mins elapsed.
## 2021-01-29 17:20:16 : Adding s1Inj to PeakMatrix for Chr (4 of 20)!, 0.143 mins elapsed.
## 2021-01-29 17:20:18 : Adding s1Inj to PeakMatrix for Chr (5 of 20)!, 0.185 mins elapsed.
## 2021-01-29 17:20:21 : Adding s1Inj to PeakMatrix for Chr (6 of 20)!, 0.226 mins elapsed.
## 2021-01-29 17:20:23 : Adding s1Inj to PeakMatrix for Chr (7 of 20)!, 0.264 mins elapsed.
## 2021-01-29 17:20:26 : Adding s1Inj to PeakMatrix for Chr (8 of 20)!, 0.305 mins elapsed.
## 2021-01-29 17:20:28 : Adding s1Inj to PeakMatrix for Chr (9 of 20)!, 0.344 mins elapsed.
## 2021-01-29 17:20:30 : Adding s1Inj to PeakMatrix for Chr (10 of 20)!, 0.381 mins elapsed.
## 2021-01-29 17:20:33 : Adding s1Inj to PeakMatrix for Chr (11 of 20)!, 0.422 mins elapsed.
## 2021-01-29 17:20:35 : Adding s1Inj to PeakMatrix for Chr (12 of 20)!, 0.464 mins elapsed.
## 2021-01-29 17:20:37 : Adding s1Inj to PeakMatrix for Chr (13 of 20)!, 0.499 mins elapsed.
## 2021-01-29 17:20:39 : Adding s1Inj to PeakMatrix for Chr (14 of 20)!, 0.535 mins elapsed.
## 2021-01-29 17:20:41 : Adding s1Inj to PeakMatrix for Chr (15 of 20)!, 0.571 mins elapsed.
## 2021-01-29 17:20:44 : Adding s1Inj to PeakMatrix for Chr (16 of 20)!, 0.606 mins elapsed.
## 2021-01-29 17:20:46 : Adding s1Inj to PeakMatrix for Chr (17 of 20)!, 0.639 mins elapsed.
## 2021-01-29 17:20:48 : Adding s1Inj to PeakMatrix for Chr (18 of 20)!, 0.675 mins elapsed.
## 2021-01-29 17:20:50 : Adding s1Inj to PeakMatrix for Chr (19 of 20)!, 0.709 mins elapsed.
## 2021-01-29 17:20:52 : Adding s1Inj to PeakMatrix for Chr (20 of 20)!, 0.745 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-01-29 17:21:01 : Adding s3Inj to PeakMatrix for Chr (1 of 20)!, 0.009 mins elapsed.
## 2021-01-29 17:21:03 : Adding s3Inj to PeakMatrix for Chr (2 of 20)!, 0.047 mins elapsed.
## 2021-01-29 17:21:05 : Adding s3Inj to PeakMatrix for Chr (3 of 20)!, 0.087 mins elapsed.
## 2021-01-29 17:21:08 : Adding s3Inj to PeakMatrix for Chr (4 of 20)!, 0.123 mins elapsed.
## 2021-01-29 17:21:10 : Adding s3Inj to PeakMatrix for Chr (5 of 20)!, 0.159 mins elapsed.
## 2021-01-29 17:21:12 : Adding s3Inj to PeakMatrix for Chr (6 of 20)!, 0.193 mins elapsed.
## 2021-01-29 17:21:14 : Adding s3Inj to PeakMatrix for Chr (7 of 20)!, 0.227 mins elapsed.
## 2021-01-29 17:21:16 : Adding s3Inj to PeakMatrix for Chr (8 of 20)!, 0.262 mins elapsed.
## 2021-01-29 17:21:18 : Adding s3Inj to PeakMatrix for Chr (9 of 20)!, 0.296 mins elapsed.
## 2021-01-29 17:21:20 : Adding s3Inj to PeakMatrix for Chr (10 of 20)!, 0.329 mins elapsed.
## 2021-01-29 17:21:22 : Adding s3Inj to PeakMatrix for Chr (11 of 20)!, 0.363 mins elapsed.
## 2021-01-29 17:21:24 : Adding s3Inj to PeakMatrix for Chr (12 of 20)!, 0.398 mins elapsed.
## 2021-01-29 17:21:26 : Adding s3Inj to PeakMatrix for Chr (13 of 20)!, 0.43 mins elapsed.
## 2021-01-29 17:21:28 : Adding s3Inj to PeakMatrix for Chr (14 of 20)!, 0.463 mins elapsed.
## 2021-01-29 17:21:30 : Adding s3Inj to PeakMatrix for Chr (15 of 20)!, 0.496 mins elapsed.
## 2021-01-29 17:21:32 : Adding s3Inj to PeakMatrix for Chr (16 of 20)!, 0.528 mins elapsed.
## 2021-01-29 17:21:34 : Adding s3Inj to PeakMatrix for Chr (17 of 20)!, 0.56 mins elapsed.
## 2021-01-29 17:21:36 : Adding s3Inj to PeakMatrix for Chr (18 of 20)!, 0.598 mins elapsed.
## 2021-01-29 17:21:38 : Adding s3Inj to PeakMatrix for Chr (19 of 20)!, 0.631 mins elapsed.
## 2021-01-29 17:21:40 : Adding s3Inj to PeakMatrix for Chr (20 of 20)!, 0.663 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-01-29 17:21:48 : Adding s3UI to PeakMatrix for Chr (1 of 20)!, 0.01 mins elapsed.
## 2021-01-29 17:21:51 : Adding s3UI to PeakMatrix for Chr (2 of 20)!, 0.046 mins elapsed.
## 2021-01-29 17:21:53 : Adding s3UI to PeakMatrix for Chr (3 of 20)!, 0.083 mins elapsed.
## 2021-01-29 17:21:55 : Adding s3UI to PeakMatrix for Chr (4 of 20)!, 0.117 mins elapsed.
## 2021-01-29 17:21:57 : Adding s3UI to PeakMatrix for Chr (5 of 20)!, 0.153 mins elapsed.
## 2021-01-29 17:21:59 : Adding s3UI to PeakMatrix for Chr (6 of 20)!, 0.19 mins elapsed.
## 2021-01-29 17:22:01 : Adding s3UI to PeakMatrix for Chr (7 of 20)!, 0.227 mins elapsed.
## 2021-01-29 17:22:03 : Adding s3UI to PeakMatrix for Chr (8 of 20)!, 0.262 mins elapsed.
## 2021-01-29 17:22:05 : Adding s3UI to PeakMatrix for Chr (9 of 20)!, 0.295 mins elapsed.
## 2021-01-29 17:22:07 : Adding s3UI to PeakMatrix for Chr (10 of 20)!, 0.328 mins elapsed.
## 2021-01-29 17:22:09 : Adding s3UI to PeakMatrix for Chr (11 of 20)!, 0.361 mins elapsed.
## 2021-01-29 17:22:12 : Adding s3UI to PeakMatrix for Chr (12 of 20)!, 0.399 mins elapsed.
## 2021-01-29 17:22:14 : Adding s3UI to PeakMatrix for Chr (13 of 20)!, 0.434 mins elapsed.
## 2021-01-29 17:22:16 : Adding s3UI to PeakMatrix for Chr (14 of 20)!, 0.468 mins elapsed.
## 2021-01-29 17:22:18 : Adding s3UI to PeakMatrix for Chr (15 of 20)!, 0.5 mins elapsed.
## 2021-01-29 17:22:20 : Adding s3UI to PeakMatrix for Chr (16 of 20)!, 0.532 mins elapsed.
## 2021-01-29 17:22:22 : Adding s3UI to PeakMatrix for Chr (17 of 20)!, 0.565 mins elapsed.
## 2021-01-29 17:22:24 : Adding s3UI to PeakMatrix for Chr (18 of 20)!, 0.599 mins elapsed.
## 2021-01-29 17:22:26 : Adding s3UI to PeakMatrix for Chr (19 of 20)!, 0.634 mins elapsed.
## 2021-01-29 17:22:28 : Adding s3UI to PeakMatrix for Chr (20 of 20)!, 0.666 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-01-29 17:22:36 : Adding s2UI to PeakMatrix for Chr (1 of 20)!, 0.009 mins elapsed.
## 2021-01-29 17:22:38 : Adding s2UI to PeakMatrix for Chr (2 of 20)!, 0.043 mins elapsed.
## 2021-01-29 17:22:40 : Adding s2UI to PeakMatrix for Chr (3 of 20)!, 0.078 mins elapsed.
## 2021-01-29 17:22:42 : Adding s2UI to PeakMatrix for Chr (4 of 20)!, 0.11 mins elapsed.
## 2021-01-29 17:22:44 : Adding s2UI to PeakMatrix for Chr (5 of 20)!, 0.144 mins elapsed.
## 2021-01-29 17:22:46 : Adding s2UI to PeakMatrix for Chr (6 of 20)!, 0.177 mins elapsed.
## 2021-01-29 17:22:48 : Adding s2UI to PeakMatrix for Chr (7 of 20)!, 0.21 mins elapsed.
## 2021-01-29 17:22:50 : Adding s2UI to PeakMatrix for Chr (8 of 20)!, 0.246 mins elapsed.
## 2021-01-29 17:22:52 : Adding s2UI to PeakMatrix for Chr (9 of 20)!, 0.28 mins elapsed.
## 2021-01-29 17:22:55 : Adding s2UI to PeakMatrix for Chr (10 of 20)!, 0.315 mins elapsed.
## 2021-01-29 17:22:57 : Adding s2UI to PeakMatrix for Chr (11 of 20)!, 0.352 mins elapsed.
## 2021-01-29 17:22:59 : Adding s2UI to PeakMatrix for Chr (12 of 20)!, 0.387 mins elapsed.
## 2021-01-29 17:23:01 : Adding s2UI to PeakMatrix for Chr (13 of 20)!, 0.42 mins elapsed.
## 2021-01-29 17:23:03 : Adding s2UI to PeakMatrix for Chr (14 of 20)!, 0.456 mins elapsed.
## 2021-01-29 17:23:05 : Adding s2UI to PeakMatrix for Chr (15 of 20)!, 0.489 mins elapsed.
## 2021-01-29 17:23:07 : Adding s2UI to PeakMatrix for Chr (16 of 20)!, 0.522 mins elapsed.
## 2021-01-29 17:23:09 : Adding s2UI to PeakMatrix for Chr (17 of 20)!, 0.553 mins elapsed.
## 2021-01-29 17:23:11 : Adding s2UI to PeakMatrix for Chr (18 of 20)!, 0.586 mins elapsed.
## 2021-01-29 17:23:13 : Adding s2UI to PeakMatrix for Chr (19 of 20)!, 0.617 mins elapsed.
## 2021-01-29 17:23:14 : Adding s2UI to PeakMatrix for Chr (20 of 20)!, 0.647 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-01-29 17:23:23 : Adding s2Inj to PeakMatrix for Chr (1 of 20)!, 0.009 mins elapsed.
## 2021-01-29 17:23:25 : Adding s2Inj to PeakMatrix for Chr (2 of 20)!, 0.042 mins elapsed.
## 2021-01-29 17:23:27 : Adding s2Inj to PeakMatrix for Chr (3 of 20)!, 0.077 mins elapsed.
## 2021-01-29 17:23:29 : Adding s2Inj to PeakMatrix for Chr (4 of 20)!, 0.11 mins elapsed.
## 2021-01-29 17:23:31 : Adding s2Inj to PeakMatrix for Chr (5 of 20)!, 0.143 mins elapsed.
## 2021-01-29 17:23:33 : Adding s2Inj to PeakMatrix for Chr (6 of 20)!, 0.175 mins elapsed.
## 2021-01-29 17:23:35 : Adding s2Inj to PeakMatrix for Chr (7 of 20)!, 0.207 mins elapsed.
## 2021-01-29 17:23:37 : Adding s2Inj to PeakMatrix for Chr (8 of 20)!, 0.241 mins elapsed.
## 2021-01-29 17:23:39 : Adding s2Inj to PeakMatrix for Chr (9 of 20)!, 0.273 mins elapsed.
## 2021-01-29 17:23:41 : Adding s2Inj to PeakMatrix for Chr (10 of 20)!, 0.305 mins elapsed.
## 2021-01-29 17:23:43 : Adding s2Inj to PeakMatrix for Chr (11 of 20)!, 0.339 mins elapsed.
## 2021-01-29 17:23:45 : Adding s2Inj to PeakMatrix for Chr (12 of 20)!, 0.372 mins elapsed.
## 2021-01-29 17:23:46 : Adding s2Inj to PeakMatrix for Chr (13 of 20)!, 0.403 mins elapsed.
## 2021-01-29 17:23:48 : Adding s2Inj to PeakMatrix for Chr (14 of 20)!, 0.435 mins elapsed.
## 2021-01-29 17:23:50 : Adding s2Inj to PeakMatrix for Chr (15 of 20)!, 0.466 mins elapsed.
## 2021-01-29 17:23:52 : Adding s2Inj to PeakMatrix for Chr (16 of 20)!, 0.498 mins elapsed.
## 2021-01-29 17:23:54 : Adding s2Inj to PeakMatrix for Chr (17 of 20)!, 0.528 mins elapsed.
## 2021-01-29 17:23:56 : Adding s2Inj to PeakMatrix for Chr (18 of 20)!, 0.561 mins elapsed.
## 2021-01-29 17:23:58 : Adding s2Inj to PeakMatrix for Chr (19 of 20)!, 0.592 mins elapsed.
## 2021-01-29 17:24:00 : Adding s2Inj to PeakMatrix for Chr (20 of 20)!, 0.622 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2021-01-29 17:24:07 : Adding s1UI to PeakMatrix for Chr (1 of 20)!, 0.009 mins elapsed.
## 2021-01-29 17:24:09 : Adding s1UI to PeakMatrix for Chr (2 of 20)!, 0.041 mins elapsed.
## 2021-01-29 17:24:11 : Adding s1UI to PeakMatrix for Chr (3 of 20)!, 0.072 mins elapsed.
## 2021-01-29 17:24:13 : Adding s1UI to PeakMatrix for Chr (4 of 20)!, 0.103 mins elapsed.
## 2021-01-29 17:24:15 : Adding s1UI to PeakMatrix for Chr (5 of 20)!, 0.134 mins elapsed.
## 2021-01-29 17:24:17 : Adding s1UI to PeakMatrix for Chr (6 of 20)!, 0.165 mins elapsed.
## 2021-01-29 17:24:19 : Adding s1UI to PeakMatrix for Chr (7 of 20)!, 0.195 mins elapsed.
## 2021-01-29 17:24:21 : Adding s1UI to PeakMatrix for Chr (8 of 20)!, 0.228 mins elapsed.
## 2021-01-29 17:24:22 : Adding s1UI to PeakMatrix for Chr (9 of 20)!, 0.259 mins elapsed.
## 2021-01-29 17:24:24 : Adding s1UI to PeakMatrix for Chr (10 of 20)!, 0.289 mins elapsed.
## 2021-01-29 17:24:26 : Adding s1UI to PeakMatrix for Chr (11 of 20)!, 0.319 mins elapsed.
## 2021-01-29 17:24:28 : Adding s1UI to PeakMatrix for Chr (12 of 20)!, 0.351 mins elapsed.
## 2021-01-29 17:24:30 : Adding s1UI to PeakMatrix for Chr (13 of 20)!, 0.382 mins elapsed.
## 2021-01-29 17:24:32 : Adding s1UI to PeakMatrix for Chr (14 of 20)!, 0.413 mins elapsed.
## 2021-01-29 17:24:33 : Adding s1UI to PeakMatrix for Chr (15 of 20)!, 0.443 mins elapsed.
## 2021-01-29 17:24:35 : Adding s1UI to PeakMatrix for Chr (16 of 20)!, 0.474 mins elapsed.
## 2021-01-29 17:24:37 : Adding s1UI to PeakMatrix for Chr (17 of 20)!, 0.504 mins elapsed.
## 2021-01-29 17:24:39 : Adding s1UI to PeakMatrix for Chr (18 of 20)!, 0.536 mins elapsed.
## 2021-01-29 17:24:41 : Adding s1UI to PeakMatrix for Chr (19 of 20)!, 0.566 mins elapsed.
## 2021-01-29 17:24:43 : Adding s1UI to PeakMatrix for Chr (20 of 20)!, 0.597 mins elapsed.
## Overriding previous entry for ReadsInPeaks
## Overriding previous entry for FRIP
## ArchR logging successful to : ArchRLogs/ArchR-addPeakMatrix-4f7ee62401e-Date-2021-01-29_Time-17-19-59.log

## get ArchR markers for motif enrichment
markerPeaks_subHybrid <- getMarkerFeatures(
    ArchRProj = archProj, 
    useMatrix = "PeakMatrix", 
    groupBy = "predAnn",
    bias = c("TSSEnrichment", "log10(nFrags)"),
    testMethod = "wilcoxon")

## ArchR logging to : ArchRLogs/ArchR-getMarkerFeatures-4f7e17b9f835-Date-2021-01-29_Time-17-24-45.log
## If there is an issue, please report to github with logFile!
## MatrixClass = Sparse.Integer.Matrix
## 2021-01-29 17:24:45 : Matching Known Biases, 0.004 mins elapsed.
## 2021-01-29 17:24:46 : Computing Pairwise Tests (1 of 5), 0.025 mins elapsed.
## 2021-01-29 17:25:12 : Computing Pairwise Tests (2 of 5), 0.458 mins elapsed.
## 2021-01-29 17:25:38 : Computing Pairwise Tests (3 of 5), 0.889 mins elapsed.
## 2021-01-29 17:26:03 : Computing Pairwise Tests (4 of 5), 1.309 mins elapsed.
## 2021-01-29 17:26:27 : Computing Pairwise Tests (5 of 5), 1.711 mins elapsed.
## ###########
## 2021-01-29 17:26:53 : Completed Pairwise Tests, 2.135 mins elapsed.
## ###########
## ArchR logging successful to : ArchRLogs/ArchR-getMarkerFeatures-4f7e17b9f835-Date-2021-01-29_Time-17-24-45.log

peakMarkers_subHybrid <- getMarkers(markerPeaks_subHybrid, 
                                    cutOff = "FDR <= 0.01 & Log2FC >= 1.25", 
                                    returnGR = TRUE)
heatmapMarkerPeaks <- plotMarkerHeatmap(
  seMarker = markerPeaks_subHybrid, 
  cutOff = "FDR <= 0.01 & Log2FC >= 1.25", 
  transpose = TRUE
)

## ArchR logging to : ArchRLogs/ArchR-plotMarkerHeatmap-4f7e32566bbe-Date-2021-01-29_Time-17-26-54.log
## If there is an issue, please report to github with logFile!
## Identified 20682 markers!
## Preparing Main Heatmap..
## ArchR logging successful to : ArchRLogs/ArchR-plotMarkerHeatmap-4f7e32566bbe-Date-2021-01-29_Time-17-26-54.log

ComplexHeatmap::draw(heatmapMarkerPeaks, heatmap_legend_side = "bot", annotation_legend_side = "bot")

Color RNA-seq data based on new markers

plotByExpr2 <- function(rna, genes, counts, name="feature", scale=TRUE){
  genes <- genes[genes %in% rownames(countHBC)]
  if(scale){
    scaledCounts <- t(scale(t(counts[genes,]+1e-5)))
    scaledCounts <- scaledCounts[!is.na(rowVars(scaledCounts)),]
    yhat <- colSums(scaledCounts)
  } else {
    yhat <- colSums(counts[genes,]+1e-5)
  }
  rna <- AddMetaData(rna, yhat, col.name = name)
  # featureplot
  p1 <- FeaturePlot(rna, features = name) + theme(legend.position = "none")
  # boxplot
  df <- data.frame(y=yhat, cluster=rna$ctHybridSub)
  p2 <- ggplot(df, aes(x=cluster, y=yhat)) +
    geom_boxplot() +
    theme_classic()
  cowplot::plot_grid(p1, p2, nrow=1)
}

peakSet <- archProj@peakSet
markersHybrid <- peakMarkers_subHybrid$C2_Hybrid
markersHybrid1 <- peakMarkers_subHybrid$hybrid1
markersHybrid2 <- peakMarkers_subHybrid$hybrid2

peaksHybrid <- peakSet[queryHits(findOverlaps(peakSet, markersHybrid, type = "equal"))]
stopifnot(length(peaksHybrid) == length(markersHybrid))

peaksHybridFunnel <- peakSet[queryHits(findOverlaps(peakSet, markersHybrid1, type = "equal"))]
stopifnot(length(peaksHybridFunnel) == length(markersHybrid1))

peaksHybridHorn <- peakSet[queryHits(findOverlaps(peakSet, markersHybrid2, type = "equal"))]
stopifnot(length(peaksHybridHorn) == length(markersHybrid2))

plotByExpr2(rna, 
           genes=intersect(unique(mcols(peaksHybrid)$nearestGene), rownames(countHBC)), 
           counts=countHBC, 
           name="HybridMarkers", 
           scale=TRUE)

plotByExpr2(rna, 
           genes=intersect(unique(mcols(peaksHybridFunnel)$nearestGene), rownames(countHBC)), 
           counts=countHBC, 
           name="Hybrid1Markers", 
           scale=TRUE)

plotByExpr2(rna, 
           genes=intersect(unique(mcols(peaksHybridHorn)$nearestGene), rownames(countHBC)), 
           counts=countHBC, 
           name="Hybrid2Markers", 
           scale=TRUE)

Motif enrichment for sub-hyrbid clusters

archProj <- addMotifAnnotations(ArchRProj = archProj, 
                                motifSet = "cisbp", 
                                name = "Motif",
                                force=TRUE)

## No methods found in package 'IRanges' for request: 'score' when loading 'TFBSTools'

## ArchR logging to : ArchRLogs/ArchR-addMotifAnnotations-4f7e1fa6f0fb-Date-2021-01-29_Time-17-27-06.log
## If there is an issue, please report to github with logFile!

## peakAnnotation name already exists! Overriding.

## 2021-01-29 17:27:07 : Gettting Motif Set, Species : Mus musculus, 0.002 mins elapsed.

## Using version 2 motifs!

## 2021-01-29 17:27:09 : Finding Motif Positions with motifmatchr!, 0.03 mins elapsed.

## 2021-01-29 17:29:53 : Creating Motif Overlap Matrix, 2.762 mins elapsed.

## 2021-01-29 17:29:55 : Finished Getting Motif Info!, 2.81 mins elapsed.

## ArchR logging successful to : ArchRLogs/ArchR-addMotifAnnotations-4f7e1fa6f0fb-Date-2021-01-29_Time-17-27-06.log

## upregulated motifs 
motifsUp <- peakAnnoEnrichment(
    seMarker = markerPeaks_subHybrid,
    ArchRProj = archProj,
    peakAnnotation = "Motif",
    cutOff = "FDR <= 0.05 & Log2FC >= 0.5"
  )

## ArchR logging to : ArchRLogs/ArchR-peakAnnoEnrichment-4f7e916fef9-Date-2021-01-29_Time-17-30-00.log
## If there is an issue, please report to github with logFile!

## 2021-01-29 17:30:04 : Computing Enrichments 1 of 5, 0.059 mins elapsed.

## 2021-01-29 17:30:04 : Computing Enrichments 2 of 5, 0.071 mins elapsed.

## 2021-01-29 17:30:05 : Computing Enrichments 3 of 5, 0.081 mins elapsed.

## 2021-01-29 17:30:06 : Computing Enrichments 4 of 5, 0.093 mins elapsed.

## 2021-01-29 17:30:06 : Computing Enrichments 5 of 5, 0.104 mins elapsed.

## ArchR logging successful to : ArchRLogs/ArchR-peakAnnoEnrichment-4f7e916fef9-Date-2021-01-29_Time-17-30-00.log

# the top 10 TFs are same for all 3 hybrid clusters which is why there's only 10 cols for all 3 clusters.
heatmapMotifsUp <- plotEnrichHeatmap(motifsUp, transpose = TRUE)

## ArchR logging to : ArchRLogs/ArchR-plotEnrichHeatmap-4f7e44bd9218-Date-2021-01-29_Time-17-30-07.log
## If there is an issue, please report to github with logFile!

## Adding Annotations..

## Preparing Main Heatmap..

ComplexHeatmap::draw(heatmapMotifsUp, heatmap_legend_side = "bot", annotation_legend_side = "bot")

# note that Gm5294 is an ortholog of the FOXL3-OT1 lncRNA and paralog of FOXL1

getTopDf <- function(motifsRes, cl){
  df <- data.frame(TF = rownames(motifsRes), mlog10Padj = assay(motifsRes)[,cl])
  df <- df[order(df$mlog10Padj, decreasing = TRUE),]
  df$rank <- seq_len(nrow(df))
  return(df)
}

plotMotifs <- function(motifsRes, cl){
  df <- getTopDf(motifsRes, cl)
  ggUp <- ggplot(df, aes(rank, mlog10Padj, color = mlog10Padj)) + 
    geom_point(size = 1) +
    ggrepel::geom_label_repel(
          data = df[rev(seq_len(30)), ], aes(x = rank, y = mlog10Padj, label = TF), 
          size = 1.5,
          nudge_x = 2,
          color = "black"
    ) + theme_ArchR() + 
    ylab("-log10(P-adj) Motif Enrichment") + 
    xlab("Rank Sorted TFs Enriched") +
    scale_color_gradientn(colors = paletteContinuous(set = "comet"))
  
  ggUp
}

# the upregulated markers for all hybrid clusters are enriched in Fox motifs.
plotMotifs(motifsUp, "C2_Hybrid") + ggtitle("Hybrid")

plotMotifs(motifsUp, "hybrid1") + ggtitle("hybrid1")

plotMotifs(motifsUp, "hybrid1") + ggtitle("hybrid1")

# Fox motifs are not enriched in either resting or activated cell types.
plotMotifs(motifsUp, "resting") + ggtitle("Resting")

plotMotifs(motifsUp, "C1_Inj") + ggtitle("Activated")

Smarcc1 motifs

Smarcc1 was recently found to regulate Keratin genes and thereby cellular fate in mammalian embryogenesis: https://www.nature.com/articles/s41586-020-2647-4. Here, we find that the markers for the activated cells are most enriched in Smarcc1 motifs. To dig a little deeper, we therefore check which of these marker peaks that may be close to Krt genes have Smarcc1 motifs here.

#detach("package:matrixStats", unload=TRUE)
# get the marker peaks for activated cells
actMarkers <- peakMarkers_subHybrid$C1_Inj

# subset to those that have Smarcc1 motif
mot <- archProj@peakAnnotation$Motif
matches <- readRDS(mot$Matches)
markerHits <- queryHits(findOverlaps(query = SummarizedExperiment::rowRanges(matches), 
                                     subject = actMarkers, 
                                     type = "equal"))
smarcHits <- which(assays(matches, "matches")[[1]][,grep(x=colnames(assays(matches, "matches")[[1]]), pattern="Smarcc1")])
markerSmarcRanges <- SummarizedExperiment::rowRanges(matches)[intersect(markerHits, smarcHits)]

# check which Krt genes are in there
krtWithMotif <- sort(unname(grep(x=mcols(markerSmarcRanges)$nearestGene, pattern="Krt", value=TRUE)))
krtWithMotif

##  [1] "Krt13"     "Krt14"     "Krt16"     "Krt18"     "Krt19"     "Krt20"    
##  [7] "Krt20"     "Krt26"     "Krt42"     "Krt5"      "Krt5"      "Krt6a"    
## [13] "Krt84"     "Krtap11-1" "Krtap2-4"  "Krtap4-16"

# plot their accessibility and expression
for(kk in 1:length(krtWithMotif)){
  gene <- krtWithMotif[kk]
  if(gene %in% rownames(rna)){
    prna <- FeaturePlot(rna, features=gene)
    patac <- FeaturePlot(atac, features=gene) + xlim(c(-6,6))
    print(cowplot::plot_grid(prna, patac, nrow=1))
  } else {
    print(FeaturePlot(atac, features=gene) + xlim(c(-6,6)))
  }
}

## Warning: Could not find Krt13 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt14 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt16 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt18 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt19 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt20 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt20 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt26 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt42 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt5 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt5 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt6a in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krt84 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krtap11-1 in the default search locations, found in
## ACTIVITY assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krtap2-4 in the default search locations, found in
## ACTIVITY assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

## Warning: Could not find Krtap4-16 in the default search locations, found in
## ACTIVITY assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

Check expression of Fox TFs and target genes in scRNA-seq data

Peaks that are more accessible in hybrid cells are upregulated in motifs bound by Fox transcription factors. We expect Fox transcription factors to be higher expressed in hybrid clusters. Since Fox transcription factors generally are repressors, we could expect the target genes of Fox transcription factors to be lower expressed in the hybrid clusters.

Evidence from the literature: - Foxg1 expression is in horizontal basal cells of OE, Foxg1 -/- mice are born with defects in the OE: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2878634/ - Fox family TF motifs are most enriched in olfactory receptors: https://arxiv.org/pdf/1201.2933.pdf - Foxl1 are markers of microvillous cells: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2878634/

# selected Fox TFs as identified with motif enrichment
# top 10 motifs TFs from hybrid cells:
selectedFoxTfs <- c("Floxl1", "Foxa1", "Foxb2", "Foxb1", "Foxc1", "Foxc2", "Foxs1", "Foxd1")
selectedFoxTfs <- intersect(selectedFoxTfs, rownames(countHBC))
countsFoxSelected <- countHBC[selectedFoxTfs,]
foxSumSelected <- colSums(countsFoxSelected)
rna <- AddMetaData(rna, foxSumSelected, col.name = "foxSumSelected")
pFoxSum <- FeaturePlot(rna, features = "foxSumSelected")

df <- data.frame(cl=rna$ctHybridSub,
                 foxSum=foxSumSelected)
pBoxFoxSum <- ggplot(df, aes(x=cl, y=foxSum)) +
  geom_boxplot() +
  theme_classic() +
  ylab("Sum of selected Fox TFs expression")

cowplot::plot_grid(pFoxSum, pBoxFoxSum)

#ggsave("../plots/foxSumRNA.png", width=8)

## scaled
scaledCountsFoxSelected <- t(scale(t(countHBC[selectedFoxTfs,]+1e-5)))
scaledCountsFoxSelected <- scaledCountsFoxSelected[!is.na(rowVars(scaledCountsFoxSelected)),]
foxSumScaledSelected <- colSums(scaledCountsFoxSelected)
rna <- AddMetaData(rna, foxSumScaledSelected, col.name = "foxSumSelected_scaled")
pFoxSumScaled <- FeaturePlot(rna, features = "foxSumSelected_scaled")

df <- data.frame(cl=rna$ctHybridSub,
                 foxScaled=foxSumScaledSelected)
pBoxFoxSumScaled <- ggplot(df, aes(x=cl, y=foxScaled)) +
  geom_boxplot() +
  theme_classic() +
  ylab("Sum of selected Fox TFs scaled expression")


cowplot::plot_grid(pFoxSumScaled, pBoxFoxSumScaled)

# ggsave("../plots/foxSumRNAScaled.png", width=8)

Composite plot

## DR of RNA 
pRNA1_fin <- pRNA1 + ggtitle("scRNA-seq data") +
  theme(legend.title = element_text(size = 12),
          legend.text = element_text(size = 10),
        legend.key.size = unit(1/2, "cm"))


## DR of RNA with ATAC markers expression highlighted
pMarkersATAC <- cowplot::plot_grid(plotlist=list(pc1,pc2,pc3), nrow=1)
pMarkersATAC

## RNA and ATAC cell type transfer
pCTrans <- cowplot::plot_grid(ppred4_2 + NoLegend(), ppred2_2, nrow=1)

## markers and motifs for subhybrid clusters
## marker
heatmapMarkerPeaksMatrix <- plotMarkerHeatmap(
  seMarker = markerPeaks_subHybrid, 
  cutOff = "FDR <= 0.01 & Log2FC >= 1.25", 
  transpose = TRUE,
  returnMat = TRUE
)

## ArchR logging to : ArchRLogs/ArchR-plotMarkerHeatmap-4f7e66eac0eb-Date-2021-01-29_Time-17-30-48.log
## If there is an issue, please report to github with logFile!

## Identified 20682 markers!

## ArchR logging successful to : ArchRLogs/ArchR-plotMarkerHeatmap-4f7e66eac0eb-Date-2021-01-29_Time-17-30-48.log

pMarker <- Heatmap(heatmapMarkerPeaksMatrix, 
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        col = paletteContinuous(set = "solarExtra", n = 100),
        show_column_names = FALSE,
        row_names_gp = gpar(fontsize = 13),
        heatmap_legend_param = list(title = "Z-score"))


## motif
heatmapMotifsMatrix <- plotEnrichHeatmap(motifsUp, transpose = TRUE, returnMatrix = TRUE)

## ArchR logging to : ArchRLogs/ArchR-plotEnrichHeatmap-4f7e3e1fb916-Date-2021-01-29_Time-17-30-49.log
## If there is an issue, please report to github with logFile!

colnames(heatmapMotifsMatrix) <- unlist(lapply(strsplit(colnames(heatmapMotifsMatrix), split="_"), "[[", 1))

pMotif <- Heatmap(heatmapMotifsMatrix, 
        cluster_rows = FALSE,
        cluster_columns = FALSE,
        col = paletteContinuous(set = "comet", n = 100),
        column_names_gp = gpar(fontsize = 11),
        row_names_gp = gpar(fontsize = 13),
        heatmap_legend_param = list(title = "-log10\n p-adj."))
pHeat <- cowplot::plot_grid(grid::grid.grabExpr(ComplexHeatmap::draw(pMarker)),
                   grid::grid.grabExpr(ComplexHeatmap::draw(pMotif)),
                   nrow = 2)
pComp1 <- cowplot::plot_grid(pRNA1_fin, pMarkersATAC, nrow=1, 
                             rel_widths = c(0.4, 1), labels = c("a", "b"))
pComp2 <- cowplot::plot_grid(pCTrans, pHeat, ncol = 1, 
                             rel_heights = c(0.8, 1), labels = c("c", "d"))
cowplot::plot_grid(pComp1, pComp2, nrow = 2,
                   rel_widths= c(1/3, 1))

# ggsave("../plots/integrationComposite.pdf", height=13, width=12)
# ggsave("../plots/integrationComposite.png", height=13, width=12)

old plots

Explore the activated HBC RNA-seq clusters

Macrophage markers

FeaturePlot(atac, features = "Msr1")

## Warning: Could not find Msr1 in the default search locations, found in ACTIVITY
## assay instead

FeaturePlot(rna, features = "Msr1")

Cell cycle genes

download.file("https://raw.githubusercontent.com/rufletch/p63-HBC-diff/master/ref/cellCycle40.txt", destfile = "~/tmp/ccGenes40.txt")
cc40Genes <- read.table("~/tmp/ccGenes40.txt", stringsAsFactors = FALSE)[,1]

countsCC40 <- countHBC[cc40Genes,]+1e-5
cc40Sum <- colSums(countsCC40)
rna <- AddMetaData(rna, cc40Sum, col.name = "cc40Sum")
FeaturePlot(rna, features = "cc40Sum")

## scaled
scaledCountsCC40 <- t(scale(t(countHBC[cc40Genes,]+1e-5)))
scaledCountsCC40 <- scaledCountsCC40[!is.na(rowVars(scaledCountsCC40)),]
cc40Sum <- colSums(scaledCountsCC40)
rna <- AddMetaData(rna, cc40Sum, col.name = "cc40Sum_scaled")
FeaturePlot(rna, features = "cc40Sum_scaled")

download.file("https://raw.githubusercontent.com/rufletch/p63-HBC-diff/master/ref/cell_cycle.txt", destfile = "~/tmp/ccGenes.txt")
ccGenes <- read.table("~/tmp/ccGenes.txt", stringsAsFactors = FALSE)[,1]
ccGenes <- ccGenes[ccGenes %in% rownames(countHBC)]

countsCC <- countHBC[ccGenes,]+1e-5
ccSum <- colSums(countsCC)
rna <- AddMetaData(rna, ccSum, col.name = "ccSum")
FeaturePlot(rna, features = "ccSum")

## scaled
scaledCountsCC <- t(scale(t(countHBC[ccGenes,]+1e-5)))
scaledCountsCC <- scaledCountsCC[!is.na(rowVars(scaledCountsCC)),]
ccSum <- colSums(scaledCountsCC)
rna <- AddMetaData(rna, ccSum, col.name = "ccSum_scaled")
FeaturePlot(rna, features = "ccSum_scaled")

Wound response genes

Genes obtained from paper on HBC injury, Fig 3D caption: `Krt6a, Krt16, Sprr1a, Sprr2a3, Krtdap, Dmkn, Sbsn, and Hbegf’

wrGenes <- c("Krt6a", "Krt16", "Sprr1a", "Sprr2a3", "Krtdap", "Dmkn", "Sbsn", "Hbegf")

countsWR <- countHBC[wrGenes,]+1e-5
wrSum <- colSums(countsWR)
rna <- AddMetaData(rna, wrSum, col.name = "wrSum")
FeaturePlot(rna, features = "wrSum")

## scaled
scaledCountsWR <- t(scale(t(countHBC[wrGenes,]+1e-5)))
scaledCountsWR <- scaledCountsWR[!is.na(rowVars(scaledCountsWR)),]
wrSum <- colSums(scaledCountsWR)
rna <- AddMetaData(rna, wrSum, col.name = "wrSum_scaled")
FeaturePlot(rna, features = "wrSum_scaled")

Respiratory markers

respGenesRebecca <- c("Gp2",    "Cyp4a12b", "Kcnj16",   "Tff2", "Chad", "Acsm1", "Kcne3", "Vtcn1",  "Kcnj15", "Cldn8",  "Pigr", "Muc5b", "Calml3", "Kng2",  "Cyp2b10")

## note that Cxcl14 and Meg3 are only described for human samples in the Methods.
respGenesDavid <- c("Foxj1", "Muc5ac", "Cxcl14", "Meg3")


## genes used by Boying
respGenesBoying <- c("Reg3g", "Muc5ac", "Muc5b", "Meg3", "Krt5", "Krt14", "Trp63")


FeaturePlot(rna, features = respGenesRebecca[1:9])

FeaturePlot(rna, features = respGenesRebecca[10:15])

FeaturePlot(rna, features = respGenesDavid)

FeaturePlot(rna, features = respGenesBoying) + NoLegend()